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PREFACE 


Durtag the past several years, rapid orbit generation techniques, based on a 
first-order application of the generalized method of averaging, have been Inves- 
tigated for the National Aeronautics and Space Administration (NASA) Goddard 
Space Flight Center (GSFC). This Investigation has culminated In the develop- 
ment of a hybrid averaged orbit generator which has been implemented In the 
Research and Development (R&D) version of the Goddard Trajectory Determin- 
ation System (GTDS). 

In order to satisfy the requirements of different audiences and because of the 
scope of this investigation, the results of the lnvestlga*'on have been dociunented 
In several parts. The primary documents are as follows: 

A Recursively Formulated First-Order Semlanalvtlc Artificial Satellite 
Theory Based on the Generalized Method of Averaging. Volume I; The 
Generalized Method of Averaging Applied to the Artificial Satellite Prob- 
lem . Computer Sciences Corporation Report No. CSC/TR-77/6010, 

Wayne D, McClain, No’'’ember 1977. 

(This document presents a discussion of the application of the generalized 
method of averaging to the artificial satellite problem; the document is 
specifically directed to the analyst.] 

A Recursively Formulated First-Order Semlanalvtlc Artificial Satellite 
Theory Based on the Generalized Method of Averaging. Volume II: The 
Explicit Development of the First-Order Averaged Equations of Motion 
for the Nonspherical Gravitational and Nonresonant Third-Body Pertur- 
bations. (The present document,. ) 

[This document presents the explicit development of the flrst-oi’der aver- 
aged equations of motion for the nonspherical gravitational and nonresonant 
third-body perturbations. The document Is directed to the analyst. ] 

System Description and User*s Guide for the GTDS R&D Averaged Orbit 
Generator . Computer Sciences Corporation Report No. CSC/SD-78/6020, 
Leo W, Early, June 1978, 

(This document presents an overview of the averaged orbit generator, a 
description of the software system, and Instructions for program execu- 
tion. The document Is directed to a general audience consisting of analysts, 
programmers, and data technicians* ] 
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JlL-xi. 


The Numerical E\'aluation of the GTDS R&D Averaged Orbit Generator . 
Computer Sciences Corporation Report No. CSC/TM-78/6138, W. D. 
McClain and L. VV. Early, September 1978 (In preparation). 

[This document is directed primarily to the anah'st and user. ) 

Status Report on Numerical Averaging Methods. Computer Sciences Cor- 
poration Report No, CSC/TM-75/6039, Anne C. Long, September 1975 

[This document presents a discussion of the numerical averaging capabil- 
ity in the GTDS R&D hybrid averaged orbit generator (parts of this docu- 
ment are superseded by the report CSC/SD-78/6020 described above). 

The document is directed primarily to the analyst. ] 

Development and Evaluation of Numerical Quadrature Procedures for Use 
in Numerically Averaged Variation-of-Parameters Orbit Generators . 
Computer Sciences Corporation Report No. CSC/TM-75/6038, Leo W. 
Early, July 1975. 

[This document is also directed primarily to the analyst. ] 

Earlier documents reporting preliminary results for both analytical and numer- 
ical averaging techniques are referenced in the above documents. 
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ABSTRACT 


This report presents, In tvvo volumes, a recursively formulated, first-order, 
semlanal^^lc artificial satellite theory, based on the generalized method of 
averaging. Volume I, which has been produced under a separate cover, dis- 
cusses the theory of the generalized method of averaging applied to the artificial 
satellite problem. 

The present volume. Volume II, presents a general first-order theory for the 
accurate computation of the long-period and secular motion of a satellite caused 
by the nonspherical gravitational field of the central body. Also, the develop- 
ment of the first-order averaged third-body disturbing function is presented, 
and the theory for the accurate computation of the long-period and secular 
motion for the special case of low-altitude nonresonant satellites is completed. 
Recursive algorithms are provided for efficient evaluation of the theory. In 
addition, several mathematical developments necessary for the construction 
of the first-order theory are presented. Also, sufficient information has been 
provided to construct the analytical formulation of the first-order short-period 
variations. 

This theory has been implemented in the Research and Development (R&D) version 
of the Goddard Trajectory Determination System (GTDS), a large orbit determi- 
nation system primarily devoted to research and development efforts supported 
by the Goddard Space Flight Center (GSFC). 
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SECTION 1 - INTRODUCTION 


The prediction and definitive determination of artificial satellite orbits Is one 
of the more computationally expensive dynamical problems today. Maintaining 
accurate ephemerides for the ever-increasing number of artificial satellites 
(which Include active scientific, defense, communication, and weather satellites 
as well as defunct satellites, launch vehicles, and other debris) requires a con- 
siderable expenditure In terms of computing time. Prelaunch mission analysis 
requires that several hundred satellite trajectories over periods of up to several 
years be generated for the purposes of lifetime and geometry constraint analysis. 

Generally, these applications fall Into two categories: those applications which 
require high accuracy, e.g, , definitive orbit determination, and the low-to- 
moderate accuracy applications referred to under the broad category of mission 
analysis. The highest accuracy requirements are obtained through the extremely 
accurate hlgh-preclslon orbit generation techniques which rely on the t xpenslve 
process of numerical integration of Newton's equations of motion or some equiv- 
alent set of differential equations. Applications with less stringent accuracy 
requirements are often treated with analytical approximations. Mission analysis 
applications are often treated with analytical approximations and. In many cases, 
even two-body mechanics Is used. 

The analytical approach to the artificial satellite problem yields a set of analyti- 
cal formulas for the coordinates or orbital elements which are usually obtained 
to first or second order in a small parameter. The approach Is to separate the 
shorr-perlod, long-period, and secular components of the motion through a series 
of canonical transformations (Reference 1). The secular contributions to the 
motion are evaluated at a given time, and the canonical transformation used to 
remove the long-period component of motion Is Inverted to provide the long- 
period motion in terms of the secular elements. Finally, the transformation 
to remove the short-period terms is Inverted and evaluated with the secular and 
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long-period contributions to the elements, thus obtaining the short-period con- 
tributions to the motion. 


Generally, existing nnah'tical satellite theories have severely restricted per- 



turbation models. In addition, these anahtlcal satellite theories are not very 
flexible with respect to the extension of the force model. This is partially 
because an analytical satellite theory has not been formulated for completely 

O 

general central-body and third-body perturbation models.*" More importantly, 
the fundamental problem of developing an accurate and flexible analytical drag 
theory remains unsolved. Consequently, as knowledge of the physical environ- 
ment (e.g. , atmospheric density, geopotential coefficients of high order and 
degree) and accuracy requirements increase, the current anah-tical theories 


rl cannot be expected to keep pace. If the more costly high-precision methods are 

to be avoided and the increased accuracy requirements are to be n.et, either 
^ , more generally formulated analytical theories must be developed, including a 

! much more satisfactory treatment of analytical drag, or an alternate approach 


I must be found. 

! 



The method of averaging offers a very promising alternative approach for the 
artificial satellite problem. This approach shares similarities with both the 
numerical high-precision and the pure analytical methods and is classified as 
a semianalytical method. In essence, this method provides the long-period 
and secular motion of the satellite very efficiently through the numerical inte- 
gration of the averaged equations of motion. In addition, the theory provides 


^Y. Hagihara (Reference 2) gives an extensive list of reierences to the work in 
artificial satellite theory. 

•> 

“Small (Reference 3) has developed a first-order analytical theory for an arbi- 
trary number of zonal harmonic terms, and Mueller (Reference 4), using the 
Poincare- Von Zelpel technique, has developed a first-order analytical theory 
for the secular and long-period motion due to an arbitrary number of zonal 
harmonic terms. Recurrence relations are used in the evaluation of both 
theories. 
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lor short-period variations in the osculating elements (Reference 5) which are I 

required for high-accuracy applications.^ 

The method of averaging approach is particularly flexible, especially with re- 
spect to the atmospheric drag models. Not only aie the complex drag models 
which are used in high-precision theories easily accommodated, but they are 
also easily interchanged without any Impact on the theory (Reference 6). Gen- 
eral models for the central body, nonspherlcal gravitational, and third-body 
perturbations can be developed in . stralghtfonvard manner, which is the sub- 
ject of this volume. 

A combination of numerical evaluation and theoretical considerations indicates 
that the method of averaging approach is generally two to three orders of mag- 
nltude more efficient than the high-precision techniques. Specifically, it has 
been shown in References 6 and 7 that a first-order application of the method 
of averagiug to the artificial satellite problem produces the long-period and g 

secular motion very accurately and with the computational efficiency cited above. | 

Thus, the method of averaging provides a low-cost, long-term orbit prediction 
capability useful for the following applications: 

• Mission analysis (lifetime and geometric constraints) 

• Tracking station acquisition schedules 

• Dynamic modeling required for differential correction (DC) pro- j! 

cedures used to solve for dynamical parameters, e.g. , high-order 1 

geopotential coefficients I 


This is equivalent to inverting the first canonclal transformation in the analyti- 
cal satellite theory to obtain first-order short-period variations In the osculat- 
ing elements, which are then superimposed on the secular and long-period 
elements. See Section 4 of Reference 5 for more details. 

9 

“For very strongly drag perturbed satellites, the Increase in efficiency may be 
reduced to between one and two orders of magnitude. 
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The effoetlveness of reprosontltts' the osculattnjj elements by superimposing first- 
order short-|)orlod variations on the mean elements has been demonstrateil by 
Lutzky and L'phoff iHeference S). Also, it can bo shown both from the tllscusslon 
In \ olume 1 of this report (Heforenee 5) and In the present volume and from the 
discussion in neference l) that the first-order short-period variations can be for- 
mulateil analytically. It appears that the cost of the evaluation of these analytical 
formulas is quite feasible and would be roughly equivalent to a single evaluation 
of the mean element rates which are numerically Integrated to obtain the long- 
period and secular motion. However, the cost of evaluating these short-period 
variations at several points In a single orbital revolution would be considerably 
less than the cost of evaluating the same number of mean element rates. ^ Con- 
sequently, It appears that this high-accuracy mode of the method of averaging 
couK! prove to be significantly more efficient than the usual hlgh-prcclslon tech- 
nhiues and, thus, may offer an efficient high-preclslon orbit prediction capability 
for tlcflnltlvc orbit determination proccilures, particularly where extended data 
Intervals with <lata g:ips are cncountereil. In addition, the short-period variations 
can ami should be useil to (.lcvelo|) an osculatlng-to-mean element conversion cap- 
ability. 


I J 


Within the same revolution or "pnss," the slowly varying mean elements are 
essentially constant. Only those functions which depend explicitly on the fast 
variable necil to be reevaluated. 



1. 1 OVERVIEW OF THE METHOD OF AVERAGING 


Tlie efficiency of the method of averaging procedure arises from the fact that 
the maximum step size which can be used in the numerical Integration of a set 
of differential equations is constrained by the highest significant frequency (l.e. , 
shortest period) contained therein. The method of averaging is used to remove 
high-frequency components from the equations of motion. The resulting aver- 
aged equations of motion are Integrated numerically but with a significantly 
gi'eater step size than can be used with the high-precision equations. The long- 
period and secular components of the satellite motion are thus obtained. The 
short-period component of the motion can be computed either numerically (Ref- 
erence S) or from aaaljTical formulas which can be constructed from the results 
contained in Reference 3 and Sections 3 and 4 of the present document. These 
formulas are also developed In Reference 9. In most cases, the computational 
savings achieved by the larger step size (which results In fewer force evalua- 
tions) far outavelghs the possible additional cost of the uerlvatlve evaluation.^ 
thereby effecting a significant decrease in the overall computational costs. 

The technique of removing the high-frequency terms from the equations of motion 
was first used by Lagrange in his investigations of the planetary motion. In the 
particular formulation of the equations of motion developed by Lagrange, the 
high-frequency terms, In the case of conservative perturbing forces, could be 
Isolated more or less by Inspection. However, a rigorous mathematical foun- 
dation for this technique was not provided until the relatively recent work by 
Krylov and Bogollubov (Reference 10) on asymptotic methods for nonlinear 
oscillations. 

Two approaches are available for the application of the method of averaging. 

The high-frequency components of the equations of motion can be removed 


The exact cost of a derivative evaluation depends on the specific perturbations 
and the characteristics of the satellite orbit, which may permit significant 
truncation of the series expansions. 





numerical^' by appllcatloa of a quadrature around an appropriate formulation 
of the hlgh-preclslon equations of motion. This procedure Is known as the 
numerical averaging approach. If the perturbing forces are conservative, the 
equations of motion can be expressed using Lagrange's formulation (Reference 5), 
and the averaging quadrature can be performed analytically. Under certain 
assumptions,^ this method produces the same result as that obtained b^' Inspec- 
tion. This semlanalytlcal procedure of numerically Integrating the analytically 
av’eraged equations of motion Is referred to as the analytical averaging approach. 


^The assumptions arise when either the Greenwich Hour Angle (l.e. , the Earth's 
rotation) or the fast variable of the disturbing third body appear In the perturba- 
tion models. Specifically, these quantities are assumed to be completely Inde- 
pendent of the satellite fast variable, both explicitly and Implicitly through the 
time. 
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1.2 RECENT DEVELOPMENTS IN ANALYTICAL AVERAGING THEORY 


Recentlj", several authors have investigated general, analytically averaged 
perturbation models for the third-body and nonspherical gravitational perturba- 
tions in terms ot nonsingular element sets. Cefola and BroucUe (Reference 11) 
developed recursively formulated models for the nonresonant third-body and 
zonal harmonic perturbations based on the nonsingular equinoctial elements. 

The development of the zonal harmonic model Is similar to that of Cook (Ref- 
erence 12), with the exception that the inclination function is developed in terms 
of associated Legendre polynomials and their derivatives and certain complex 
polynomials. Cefola's third-body model is developed in terms of the direction 
cosines of the disturbing third-body position vector, which proves computation- 
ally efficient but is limited to nonresonunt cases, Cefola outlined an extension 
of his zonal harmonic model to Include the nonresonant tesseral harmonic terms 
(Reference 13) and later completed and extended the model to include resonant 
phenomena (Reference 14). Giacaglia (Reference 15) reformulated Kaula's 
(References 16 and 17) perturbation’ models (using Allan's inclination function) 
in a nonsingular element set and provided a set of recursive algorithms for 
computational purposes. Finally, Nacozy and Dallas (Reference 18) also 
reformulated the Kaula geopotential model (using Allan's inclination function) 
in terms of a nonsingular element set. No recursive al(tovithms were provided. 




1.3 SURBURY 

This report is the result of a series of task assignments with the objective of 
Imr' ' mentlng In the Research and Development (R&D) version of the Goddard 
Trajectory Determination System (GTDS) a set of recursively evaluated, first- 
order analytically averaged equations of motion for an artificial satellite per- 
turbed by nonresonant third-body and nonspherlcal gravitational perturbations. 
This analtylcal-averaglng capability enhances the GTDS numerical averaging 
capability (Reference 6) and provides for optimal averaged perturbation models 
for each specific type of perturbation (with the exception of third-body resonance 
cases, which were not considered). Partial results obtained for some of the 
otpimal averaged perturbation models In GTDS have been presented In Refer- 
ence 7. 

For Implementation, Cefola’s averaged perturbation models (Reference 11) for 
the nonresonant third-body and zonal harmonic perturbations are adopted. The 
nonresonant tesseral harmonic model was developed as part of this task assign- 
ment using the approach outlined by Cefola in Reference 13. The resonant 
tesseral harmonics model was also developed as part of the task assignment 
from a completely general nonspherlcal gravitational theory designed to yield 
the zonal harmonics, nonresonant tesseral harmonics, and resonant tesseral 
harmonics models as special cases. In addition, all models were generalized 
to handle retrograde as well as direct equinoctial elements (see Appendix A of 



Reference 5). 

The brute-force Implementation of recursive algorithms can contribute to com- 
putational Inefficiency and can possibly Introduce artificial singularities (not In 
the equations of motion, but In the model evaluation). To Insure against this 
possibility, careful consideration was given to the ordering of the terms In the 
models, such that the recurrence formulas proceed In the proper direction to 
avoid small divisors. Also, the amount of recomputation and storage require- 
ments are minimized. 
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For Implementation In the GTDS RiD system, It was felt that the resonant tes- 
seral harmonic motJel should be very flexible with respect to the specific reso- 
nant harmonic terms used. Tlie existence of a resonance dictates which terms 
In the potential expansion are significant to the loug-{)erlod motion. Knowledge 
of the common characteristics of these terms and the proper use of the recur- 
sive algorithms could have provided a means for further optimization of this 
model. However, the procedure would have been automatic, with the program 
expecting a certain set of terms. Therefore, for the purposes of flexibility and 
at some additional comp, tatlonal costs, the contributions from each spherical 
harmonic term are computed entirely Independent^' of all other terms. ^ 

Due to the extensive new softwai'e for the analytical averaging capability, as 
well as to the extensive modifications required to the prevlousl^v Implemented 
averaging software (particularly’ the Input processor and Initialization procedures 
and the attendant added complexity of executing the GTDS R&D averaging capabil- 
ity), It was decided that a system description and user's guide for the GTDS 
R&D averaging capability would bo issued under a separate cover (Reference 19), 
In addition, a document extending the numerical results beyond those pi'esented 
In Reference 7 Is also In preparation. This document (Reference 20) will 
discuss the computational costs In terms of machine processing time, the ac- 
curacy of the analytical averaging methods, and the procedure and algorithms 
used to develop an automatic truncation capability to further optimize the per- 
turbation models for each particular case. 

The present report consists of two volumes. \'olume 1 (Reference 5) presents 
a comprehensive discussion of the application of the generalized method of aver- 
aging to the artificial satellite problem and the resulting formulation of the 
averaged equations of motion. Included In the discussion are the formulation 
of the N’arlatlon of Parameters (VOP) equations of motion and the application 


^The capability to automatically' select the resonant terms was Implemented In the 
GTDS R&D version. However, no special relationship among them Is assumed. 
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of the method of averaging to the \’OP equations of motion. Other topics dis- 
cussed Include the criteria for the selection of short-period terms, the applica- 
tion of the method of averaging to the case of two or more perturbing functions, 
the application of the method of averaging to cases involving resonance phenom- 
ena, and a discussion of the first-order short-period variations in the osculating 
elements and their application to osculating-to-mean and mean-to-osculating 
element conversions. 

\'’olume II (the present document) presents the mathematical formulation, in 
nonslngular equinoctial elements, of the nonsphertcal gravitational and non- 
resonant third-body disturbing functions required for the first-order averaged 
equations of motion. Section 2 of this document presents some mathematical 
developments required for the expansion of the disturbing functions. Specifically, 
Section 2. 1 discusses the theory of the rotation of spherical harmonic functions. 
Next, Section 2, 2 develops certain Fourier series expansions which are of im- 
portance in the development of the disturbing functions. 

Section 3 presents the explicit theor}' [or the nonspherical gravitational pertur- 
bation. The development of the nonspherical gravitational disturbing function 
is discussed in Section 3. 1, and the disturbing function is expressed in equinoc- 
tial elements In Section 3.2. Also, a discussion relating Kaula's Inclination 
function (References 16 and 17) to the inclination function developed In this 
report is presented. The nonspherical gravitational disturbing function is aver- 
aged In Section 3.3. The averaging operation and the concepts and implications 
of time-dependent and time-independent averaging are discussed. In addition, 
the averaged disturbing functions for the special cases of the zonal harmonics, 
combined zonal and nonresonant tesseral harmonics, and resonant tesseral 
harmonics are Isolated and presented. In Section 3. 4, the partial derivatives 
of the nonsphertcal gravitational disturbing function which are required for the 
av'eraged equations of motion are presented for each case, and the recurrence 
relations used for the evaluation of the constituent functions are given. 
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Sectloa 4 of the present volume presents the exi'Ueit theory for the disturbing 
third-body perturbation. Section 4. 1 discusses the development of the third- 
body disturbing function, and Section 4. 2 gives the general expansion of the 
disturbing function. Due to time and other resource constraints, onl>' an outline 
of the general development is presented. However, because of the similarities 
with the nonspherical gravitational theory presented in Section 3, the neglected 
details are straightforward. Section 4.3 presents another expansion of the 
third-body disturbing function which is well suited for cases of nonresonant 
(with the third body) near-Earth satellites. Section 4.4 presents the partial 
derivatives of the averaged disturbing function developed in Section 4. 3 which 
are required for the averaged equations of motion in the special case. The 
necessary recurrence relations for evaluation of the theory are also provided. 

The equations of motion for all models arc given in what is considered to be an 
optimal form, taking into account the minimization of the combined computa- 
tlona! and storage costs while avoiding computational singularities. It is this 
final form of the models tliat was implemented in the GTDS R&D system. These 
models reflect, to some extent, the computer environment in which they were 
implemented, i.e. , the GSFC IBM S, 360-95 computer. A variable can assume 
a moderatel>‘ wide range of magnitudes in this environment. This, of course, 
is not true in all computer environments. Thus, for Implementation on other 
computers, it may be necessary to Introiiuce normalization factors and to other- 
wise redefine certain functions appearing in the models given in Sections 3 and 
4 in order to minimize the magnitude differences between quantities. 
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1.4 CONCLL'SIONS 

Most satolllto theories bused on an avoraglnif concept, e.g. , Cefola's (Refer- 
ence 14), GlacagUa’s (Reference 15), and Kaula's (References Iti and 17), have 
been formulated using the classical assumptions of a time-independent pertur- 
bation model^ or of exact resonance. While these classical assumptions guar- 
antee ■'**» removal of tlie unwanted short-period contributions to the motion, 
they mai‘, in some cases, produce a significant exaggeration of certain medlum- 
r.nd long-period contributions to the motion. TlUs fact should be considered in 
the implementation of an averaging theory. 

The zonal, nonresonant tesseral, and resonant tesseral harmonic models, 
which are presented in Section 3 of this volume and which have been implemented 
in the GTDS R&D system, comprise a completelj' general first-order theory for 
the contributions to the long-period and secular motion caused by the central- 
body nonspherlcal gravitational field. However, as discussed in Section 3. 1 of 
Reference 5, it is recommended that the effects contributed by the nonresonnnt 
tesseral harmonics be excluded from tlw averaged equations of motion, since 
they unnecessarll.v restrict the integration step size and since they can be for- 
mulated analj ticalVv in the same way as the first-order short-period variations 
in the osculating elements. In addition, the analy tical formulas for the first- 
order short-period variations in the osculating elements can. In essence, be 
developed from the information contained in both volumes of this report. In 
order to obtain the* final formulas. It is only necessary to extract the appropri- 
ate formulas from Sections 2* 2, 3, and 4 of the present v'olumc and substitute 
them into Equations (4-14) of Section 4 in \'olun\e 1. 

Numerical evaluation of the nonspherlcal gravitational theories for long-period 
and secular motion, |)erformed as part of the investigation and documented in 


See Section 3.3 for more details. 
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Hefereaces 7 and 20, indicates that the theory is very efficient and accurate for 
all cases considered, except resonant effects on larjfe-eccentrlcity orbits. The 
pro|)er treatment of the eccentricity expansions (Hansen coefficients) in these 
cases still remains an open question. For the present, the numerical avcrairing 
method, with a properl>’ reduced force model, provides au acceptable alternative 
for cases of deep resonance* However, the effects of the high eccentricities are 
observed through higher order quadrature algorithms. 

Although the general third-body disturbing function is developed in Section 4.2 
of the present volume, onlj- the special case discussed in Section 4. 3 was im- 
plemented in the GTDS R&D s.vstem. This implementation is restricted to those 
applications involving low to moderate altitude satellites that are not in reso- 
nance with the disturbing third body. For Eartli satellites, these requirements 
ti'anslate into satellites with periods less than 3 to 4 days and which are not in 
resonance with the ^loon. 

The numerical evaluation of this capability (References 7 and 20) indicates that 
it produces the long-period and secular motion of nonresonanf satellites of low 
to moiierate altitude very accurateli' and efflclentl\'. This theory is comuletel.v 
inadequate for Earth satellites with periods longer than 4 days. Also, the 
expense of numerteaH,v averaging such cases is completel.v unacceptable since 
• it is at least as expensive, and usuall>- more expensive, than high-precision 
techniques. 

The application of averaging methods to satellites with longer periods requires 
a double-averaged (nonresonant cases) third-bod.v theory or a single-averaged 
third-body resonance theory, depending on the orbital characteristics of the 
satellite. It also seems feasible for such theories that first-ord?r short-period 
variations in the osculating elements can be anuU'ticall.v formulated to meet 
high-precision requirements. 
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SECTION 2 - MATHE^IATICAL PRELIMINARIES 


The explicit development in equinoctial elements of the nonspherical gravitational 
and third-body disturbing functions requires the mathematical formulation for the 
rotation of the spherical harmonic functions. In addition, Fourier series repre- 
sentations are required for functions of the form 

(L^ ( COS sL I 
W) I iin 

where r is the radial distance of the satellite, d is the semimajor axis of the 
osculating orbit, L is the true longitude of the satellite which describes the sat- 
ellite position in the orbit relative to the origin of the longitudes, and s is an 
integer. 

This section presents a general discussion of the rotation of the spherical har- 
monic functions and develops the Fourier expansions, in the true, eccentric, and 
mean longitudes, of the functions 



where 

jsL 

C a dxp(j%L) • COS sL ♦* ] Sin sL 

The results of this section are applied in Sections 3 and 4 to obtain the disturbing 
functions and averaged equations of motion for the nonspherical gravitational and 
third-body perturbations, respectively. 


2. 1 ROTATION OF THE SPHERICAL HARMONIC FUNCTIONS 


A spherical harmonic function takes the form 


(-M 


sm m. 


where (r, A»^) are the spherical coordinates of a given point, i.e. , the satellite 
position, and is the associated Legendre polynomial of degree X and 

order m and is defined to be (see Section 3. 1 for more details) 


*n/i 1 a vX 

^ MSXSI) (2-2) 


The quantities Cj^^ and S^ ^ are referred to as the spherical harmonic coef- 
ficients and are usually empirically determined constants, 

A complex variable representation of the spherical harmonic function is useful for 
the purposes of this discussion and Equation (2-1) is expressed in the form 


I 

'** I -M ( • i Km) Pi, 


where Re designates the real part and J is the '.maginary unit, i. e. , 


j • 


Often, it is necessary to transform the above expression to a different reference 
system in order to obtain an expression in terms of the transformed coordinates. 
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Most generally, such a transformation involves a t.> anslation and rotation of the 
coordl'iate reference system.^ If, however, the two coordinate rcfcieacc systems 
possess a common origin, the transformation reduces to only a rotation and the 
spherical harmonic function is to be expressed in the transformed latitude and 
longitude (0*, A'). Since the radial distance, r, is invariant under a rotation 
and since the spherical harmonic coefficients are independent of the position, it 
is sufficient to rotate onl>* the surface harmonic function 

enp (jmA) (2-1) 

In this section, a discussion of the mathematical formulatior ^ a general rotation 
is presented. Specificall.v, the Euler angle and Euler parameter representations 
of a rotation are presented. Next, the rotation of a surface harmonic Is discussed. 
I^rticular attention is devoted to the development of the generalized inclination 
function in terms of orthogonal pol> nomials. 

2. 1. 1 Mathematical Representation of a General Rotation 

Since a general rotation of a coordinate reference system leaves the origin invar- 
iant, onl,v three independent parameters are required to describe the rotation, 
tt follows that these three independent parameters specify the relative orientation 
of two coordinate reference systems with a common origin. The three parame- 
ters most frequentl>’ cHosen are the Euler angles 0Sn<2w, 0StJ12ir, and 
OSi^tr shown in Figure 1. 

2. 1. 1. 1 The Euler Angle Representation 

The primed coordinate system in Figure 1 cun be obtained by performing three 
simple rotations of the unprtmed system. Specifically’, the first rotation is per- 
formed about the unprimed z axis through the angle A * The second rotation is 

^Lee (Reference 21> discusses the effects on the spherical harmonic furctions 
caused by this more general transformation. Also, see Aardoom (Reference 22). 
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Figure 2-1. Euler Angles 



performed about the new x axis through the angle I , thus rotating the xy plane 
into the x*y' plane. Flmlly, the third rotation is performed about the z’ axis 
through the angle OJ . The general rotation is then expressed as a composite of 
the three simple rotations, l.e. , 


where 


and 


A « R^(i) R^ia) 


t 0 

« 0 COS0 

0 - i\r\ & 

COS 9 sm 9 0 

R^(9) * - sin 9 COS 9 0 

0 0 1 



(2-5) 


(2-6) 


(2-7) 
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are the matrix representations of rotations about the x and z axes, respectively. 
The product of the three rotation matrices in Equation (i-5) yields the following 
general rotation matrix: 




Sy - Cj 


“ ^w^rx 

- Si Co. 


C^,S; 


where - 

» C0«» X 

• s'm 

This transformation matrix is used to transform the coordinates of the position 
vector in the unprimed reference system to the coordinates in the primed refer- 
ence system through the equation 


r' . 


Tlie inverse transformation 


f - A'^ 


is easily obtained since the transformation is orthogonal, i.e. , 


Z ^\\ 


(k=l,2,3) 


where the quantities m elements of the transformation matrix A , i, e. , 


A ' I n„,, 


(2-12) 
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Maktng the substitutions 


O* a 


Cl" Q 

\ 


(2-16a) 


CImQ 

P ' ~ 


(2-lGb) 


I t ' 

I 1; 


It follows that 


where 


A*^ * A'^ 


A * U 


I “*rv^,n 


(2-13) 


(2-14) 


is the ti'aiispose of the A matrix obtained b\’ Interchanging the rows and columns 
of the matrix (Reference 23). 

2. 1. 1. 2 Tlie Euler Parameter Representation 

Other representations of the transformation matrix exist. For example, by 
using elementary trigonometric identities. It Is easily shown that matrix A can 
be expressed as 



^'t/A 


Su,Si 


Si 

a. 
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and using the identities 


(2- 16c) 
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L' 


Ca. • - S‘ 


(2-17a) 
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1 yields the matrix expression 
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are the well-known Euler parameters (Rofei'ence 24). 

llie Euler parameters are frequently used in the quaternion representation for 
the transformed coox'dlnates, i.e. , 


r' . (J? H**- 


(2-20) 


where the quaternion q is a hj'percomplex variable defined to be 


'.nd where 


It' *vi "Ij'' 


t » XI ^ j 4* ik 




• -1 


i ] * - j i « k 


jk « -k^ s i 


k\ » -\K a j 


Equation (2-20) should be taken onl,v as a convenient algorithm and does not imply 
that the quaternion is the rotational operator (Reference 24). 
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2.1.2 Rotation of a Sm’face Harmonic Function 


Essentially, thei'e are two approaches to the rotation of the surface harmonic func- 
tion 





(2-21) 


One approach relies on the brute-force substitution of the transformed coordinates 
using the expressions^ 


4in (P « 


1 - f 



e 



Cj 


( 2 - 22 ) 
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Wa. * 



-a 



S; e 


j(H' ir/a) 




(2-23) 


Into Equations (2-21). Equations 02-22) and (2-23) are obtained from Equations 
(2-10), (2-13), and . (2-8) or (2-151. This approach can yield complicated expres- 
sions which may obscure the real nature of the ti'ansformatlon. 

The second approach, based on the work of G. Herglotz and W, Magnus (Refer- 
ence 25) Is to rotate the entire expression for the surface harmonic instead of 
direct substitution of the rotated position components. Since the derivation of 
the theory Is quite lengthy and since It Is provided In Reference 25 and also by 
Lee (Reference 21), it will not be presented here. 

^In real variables, Equation f2-22) takes on the more familiar form 

sin ^ « sin i sin ( A* ♦to) cos ♦ cosi sin<;^i' 

The real and Imaginary parts of Equation (2-23) also assume a familiar form. 
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The rotation of a surface harmonic function of (A, <<>) to a reference system with 
the corresponding coordinates 4>') takes the form 




Sa-i. 


(2-24) 


where the associated Legendre polynomial of negative order is defined in terms 
of the corresponding polynomial of posltlvie order by the relation 


CO 

X,-n 



(L-Ol 

71 — T7 h «Cx^ 
(JL+nV. 


(n>0) (2-25) 


ni s 

The function {p, <r, T) Is discussed next, 

2. 1.2, 1 Development of the Function ^(p, <T,X) 
ni s 

Tlie function (p,a,Z) is defined to bo 


m,s { r / -TT \ 1) 

SjiCyO,cr,T) * e^pp (m-s)|op- -j - j U^^(t) (2-26) 


where 

r 

(f;:) 

1 Cx 5^ Cx} 

(for m+s SOi 

(2-27a) 

- < 

1 

1 

1 

li-*; 

k ltH% %*m / 3 y 

) Cx ^ Cxy 

(for m+s > 0) 

(2-27b) 


and corresponds, in part, to the inclination function In References 16 and IT. 
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The parameters /O, (T.r describe the orlentatloa of the unprtmed or original 


coordinate reference system relative to the primed or transformed reference 
system . The corresponding transformation from the primed reference system 
to the unprlmed reference system Is given by Equation (2-9) If the symbols "r 
and f' are Interchanged. The parameters yO « C, and T are defined In terms 
of the Euler angles through Equation (2-16).^ Otherwise, the transformation 
given In Equation (2-10) Is used and the Euler angles of the inverse transforma- 
tion (fl*, I*, (J') are required to determine yO, <T, and Z through Equations (2-16). 

The notation F(a, b, c, x) In Equations (2-27) designates the well-known hyper- 
geometric series defined by (Reference 26) 




(c)n 


where Pochhammer’s symbol, (ft)jj, Is defined by 

(o.)„ * a.(oL^l)(a.'*‘a)”‘ (cL+n-i) 

Clearly, If <L Is a nonnegative Integer, then Equation (2-29) can be expressed 
using the well-known Gamma Function as 

‘ " r(o.) 


The resulting definition of o* differs from that found In Courant and Hilbert 
(Reference 25) by 7T/2 , I, e. , 

*rr 

In addition. Equation (2-26) Is a modification of the corresponding equation In 
Reference 25 In order to account for the change In the definition of a* . 
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For negative integers a, Equation (2-29) can be expressed as 


/ X / . \n C* a) • 

U\ = C-0 ; 

^ (-a-r\V 


(2-31) 


Ultimately, a set of recursive algorithms is desired for evaluating the function 
tn s 

{p, or, z) . Gauss' contiguous relations for the hypergeometric series 
(Reference 26) could be used as a basis for these recurrence relations; however, 
they are not well suited for the purposes of this investigation. This is discussed 
in Section 2. 2. 1. 4. 

Inspection of Equations (2-28) and (2-29) indicates that if either of the first two 
arguments is a negative integer, the hypergeometric series terminates in a 
polynomial. These polynomials are the orthogonal polynomials named after 
Jacobi, and they possess some very useful recurrence relations (as will be 
shown in Section 3). 

m s 

2. 1.2.2 Jacobi Polynomial Representation of the Function S^j^’ (p,<T,Z) 

Inspection of Equations (2-27) indicates that the hypergeometric series termi- 
nates to yield a polynomial of degree Is | . The relationship between the 
hypergeometric series and the Jacobi polynomial taken from Reference 26 Is 


F(-w, x) ■ 


(2-32) 


where the Jacobi polynomial takes the form 




(2-33) 
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or, equivalently, 


p».b( , ru>n^i) — y ,„, r_u.h««.m>i) , 

" niru.wi) Z^'"' rru.m.q 


(2-34) 


In principle, the Indexes a, b can assume any value with the exception of nega- 
tive Integers, l.e. , 


d^~n i b#-m 


(2-35) 


Howexer, for this Investigation, only nonnegative Integer values are considered* 

Applying Equation (2-32) to the hypergeometrlc series In Equations (2-27) and 
noting that 

i- ld\ • - Cl, 


yields 


F(-i-5, ci) . pr‘'"'%c. ) r- 

’ U-m)l ^ ' 


for m + s :S 0 and 


. - (1117^^' P. Car) (2-37) 


for *n + 8 > 0 . 
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Care must be taken to insure that the Jacobi polynomial is valid over the entire 
range of validity for the hypergeometric series. Since only integer values of m 
and s are of concern. Equations (2-35) can be expressed as 

CLiOi b iO 

« 

For Equation 0^-36), the constraints on the Jacobi polynomial are 

-nn-s i 0 5 m-s > 0 

which are simultaneously satisfied only by those values of s where s < -m . 

Hence, the Jacobi polynomial and hypergeometric series in Equation (2-36) are 
valid over the same range, i. e. , m +s 0 . The constraints on the Jacobi poly- 
nomial in Equation (2-37) are 

i 0 } s-m hO 

which are simultaneously satisfied only for sim. Thus, while the hypergeo- 
metrlc series is valid for m + s >0 or s> -m, the Jacobi polynomial is valid 
only for s > tn . 

A valid Jacobi polynomial representation for Equation (2-37) over the range 
-m ^ a 1 m is obtained through a linear transformation of the hypergeometric 
series (Reference 26), i. e. , 

li*ryi+s i ^ C^) (2-38) 

Both of these hypergeometric series are valid when m -t-s 20 is satisfied. It 
follows from the definition of the hypergeometric series that the first two argu- 
ments can be interchanged as follows: 

m-A, C-p ) * F(m-£, j C-c) (2-39) 
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In view of Equation (2-32), 
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The constraints on the Jacobi polj’nomial are m + s >0 and m - s >0 , which are 
satisfied simultaneously by -m < ? ^ m . In summary, the hypergeometric series 
in Equation (2-27) can be expressed as 


F(»*l, j C^) 


(l-w)! U,**!! P, (-Cj,) (-m i s ^ m) (2-4 la) 

Cits')! 




mtS . 

(-Ca») (mss^l) (2-41b) 


Substituting Equations (2-36) and (2-41) into Equations (2-27) and using the rela- 
tion 
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(mSs<^) (2-43C) 
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Substitution of this expression into Equation (2-26) yields the complete expres- 

in s 

Sion for the function 0’,r) In terms of the Jacobi polynomials, l.e. , 




(-0 S, P„. (C»,1 


, , (l-»w)! «’*♦**"*•*«•*•'•*♦•,* ^ 


G-r P,., (CiO 


(-i£si-m) (2-44a) 

(-m i s i m) (2 -44b) 
(m£sSi) (2-44c) 


m s 

2. 1,2.3 The Function in Terms of the Euler Angles 

Expressing Equations (2-44) explicitly In the Euler angles (fl, ui, 1), which describe 
the I'otation from the primed reference system to the unprlmed reference system, 
through Equations (2- 16) yields 
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(-1^8^-m) (2-43a) 

(-m S 8 S m) (2-45b) 


(m i s S 1 ) (2-45c) 


2-16 





4 -^4 --4 v] -4 .4 j f . r 


I ■ t 


- i 


I 

I 

e 

0 


?! i : 


|Q 

il 

i 


\ ■ 


i > 




M I ,= 

■ ■ y 

i ' ^:l 




I i 


11 

a 

a 

u 


: I n 

I a 




Lj 

u 

u 

I 

I 

I 


L.U. 


5 » 


1 t \ y ' ^ 1 t " ^ ^ f 

’. f ; ■■ ;■?> ■' ■■ ! *’ -. 




ni s 

:2, 1. J.4 Symmetry of the Funetlon * 


OK’U’INAT, PAGR IS 
Oh' I'u'*' QU.VUTY 


In view of the symmetry In Equation (--• lo), i* ie- possible to formally collapse the 
definition of the funetlon. Defining 6 = 1 for s <0 (4 - *1 for s>0) and 
using the relation given in Equation 02-42) yields the following erjpresslons for 
Equations (2-4oa) and (2-45c); 


-j(w-4»)x |«.m «v*m »*m,nw» 

* ‘ e • Cirt Si, I Pj., UC-) 


( 2 - 46 ) 


for 4 = - 1 and s > m . 

Similarly, Equation (2-45b) is symmetric about s = 0 and can be exTJressed as 


e * € (*i) € 


for 4 = - 1 and 0 S s S m . 

Tlie rotation of a surface harmonic function given by Eqi-..;.ion 0-"-4) can then be 
expressed as 




Jb 

$m0 


issA* 


(2-4S) 


. here 4 takes on the values 1 1 , 


K • I 


I 1/a. 4-0 
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s ^Mx 


(m<s^i) (2-50b) 


In view of Equation (2-25) 
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£ rr PpstO 

(ItsV. 


(s > 0) (2-51) 


Thei'efore, Equation (2-48) can be simplified to read 
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(2-52) 


for 6 = * 1 » since 
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* (5-i')l 


(2-53) 


Equations (2-50) and (2-52) are the recommended formulation for those cases where 
the symmetry is easily taken advantage of. Otherwise, the formulation should be 
developed using Equations (2-24) and (2-45). 
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2.2 EXPANSIONS OF THE PRODUCT (r/o.)" 
Expansions of the form 



(2-54) 


(where r, a, and L are the radial distance, semlmajor axis, and true longi- 
tude, respectively, and X Is the true (L), eccentric (F), or mean (A) longitude) 
play a major role In the development in equinoctial elements of the disturbing 
functions for the nonspherical gravitational and third-body perturbations. Con- 
sequently, they are also Important in the development of the analj'tlcally averaged 
equations of motion and In the analytical development of the short-period varia- 
tions In the osculating elements. 

For certain cases, each longitude possesses a particular advantage as the expan- 
sion variable, X. Specifically, for n< 0 , the above expansion Is finite In terms 
of the true longitude, and, for n>0 , the expansion Is finite In terms of the eccen- 
tric longitude. Wlille the expansion In the mean longitude Is always Infinite, It is 
of considerable Importance because of the simple relationship between the time 
and mean longitude. 

Similar expansions of the form 



V 


(2-55) 


(where f designates the true anomaly and x is the true (f), the eccentric (u), or 
the mean (i) anomaly) played an Important role in the development of the classi- 
cal disturbing function of planetary theory. These expansions were investigated 
extensively by Hansen (Reference 27) and good, if less exhaustive, discussion 
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can be found In Brown and Shook (Reference 28). A somewhat more theoretical 
discussion, based on Cauchy's first and second theorems, is given by Haglhara 
(Reference 29). 

2.2.1 Reduction to the Expansion of the Product (r/Q,)° 

Hansen's results are directly applicable to the expansions In the longitude (Equa- 
tion (2-54)). This Is easily demonstrated using the relation between the equinoc- 
tial longitudes and the corresponding classical anomalies 

X * + Ci3 in . (2-56) 

In view of Equation (2-56), the left-hand side of Equation (2-54) can be expressed 
as 


/jr)” ]sL / js (u5-»-Ifl) jsf 

\(l' ^ * 111) ^ ^ (2-57) 


Substituting the expansion In Equation (2-55) yields 


r jsL js(oJ4in) 
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Z n,s jt: 

a, e' 


(2-58) 


Inverting Equation (2-56) and substituting the result Into the expansion in Equa- 
tion (2-58) yields 


/r\'’ ]%l jsCw+in') ^ n,s It(X-u)-lfl) 


I 




(2-59) 


i:X 


I 

fi 

0 

II 
0 
c 

li 

I! 

11 

a 

II 

0 

a 

0 

0 


ii 


2-20 



i-_.. 




-f 
i-4 


4 *^'1 


I 


■i 

I 


• r : 

\ 

) 

mkm 


t 




v ; 


A comparlsou of Equations (2-54) and 02-59) indicates that 


j(^*t^(u)v iri) « t 

A ^ = e OL ^' ( 2 - 60 ) 


In view of the definition of the equinoctial elements h, k (Reference 5, Appendix A), 


and 


j (u) + in) , 

e « & ( u 1- jVi ) 

-](u)»rn) .1 

6 • e (k-]h) 


(2-6 la) 
(2 -6 lb) 


Therefore, 


or 


' i’UiK' Al.irv 


A^*€. (k-t-jli) <x^ 


,n,s» ■ ■S't . 

“ e (k-^U) 


(2-62a) 

C2-62b) 


ami the applicability of Hansen's results is demonstrated. 


The explicit development of the Fourier series expansions of the form .vflven In 
Equation (2-55) Is presented next. Much of the discussion follows the approach 
of Hansen but Is of much more limited scoiw. In addition, some results obtained 
by Hill and Newcomb, for expansions In tlie mean anomal^v, will be presented. 

2,2, 1, 1 Expansion In the IVue Anomal.v 


The Fourier expansion of the form 



(2-63) 
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is desired. I\Iore correctly, the coefficients in the expansion have only two 
indexes since the Fourier expansion of the imaginary exponential function is 


not required. If 


(fr*Zv:ei' 


(2-64) 


/ r\’^ isf r 

(t) ^ V, 


n iCl+sH 

. e 


(2-65) 


and it is sufficient to develop the expansion for (r/a)“ (Equation (2-64)). 
It follows from Fourier’s Theorem (Reference 30) that 


V. = 


it 

a. / 


(2-66) 


This expression is obtained by multiplying 6 ^ through Equation (2-64) and 
integi’atlng both sides appropriately. The resulting series collapses to a single 
term, specifically, the term when t = k, in view of the orthogonality conditions 
for Fourier series or, equivalently, the 27T periodicity of the imaginary expo- 
nential function. 

Since r/CL is a real function and since Equation (2-64) is also satisfied by the 
complex conjugates of each side, it follows that 
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A compai'lson of Equations (2-64) and (2-67) shows that 


V a V 

-t 


(2-68) 


An explicit representation of the coefficient can be obtained by a brute-force 
expansion of (r/a)*^ using the well-known relation 


1- e 


OK <2-69) 




Using Hansen's approach, the following definitions are made: 


SiVi * e 


(2-70) 


X * e 


(2-71) 


Then, 


COS a V 1 - e’ 


(2-72) 


^ a ian (<^/a) 


1 + Vi-e*- 


(2-73) 


and Equation (2-69) can be expressed as 









t I 



Therefore, 

('^) * (l+/3Xr(l+ (2-75) 

{^) * (1-/3^) tiS (i+/3X) y ) (2-76) 

for Q, a noanegatlve Integer. Hence, the expansions for (r/a)*“ reduces to the 
expansions of the product 



A/ 


(2-77) 


2. 2. 1. 1. 1 Expansions of the product (1 ^/Sx)”*^ (1 + (|3/x)]”*^ 
Using the Binomial Theorem yields the following result; 
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(2-78) 


If the following definitions are made: 

t » k-m 
p » k+m 


(2-7 9a) 
(2-79b) 
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It follows that -00 < t < 00 , since 0 $ k < oo and 0 < m < oo . To determine the 
range of p. It Is necessary’ to Invert Equations (2-79) to yield 


1 Li 


pi-t 

k . — iO 


Clearly, pi 1 1 | and p±t must be even; consequently, p is defined as 


p » -h |t 


Then, if 1 1 ( £ p < oo , it follows that 0 < I < oo and 


rvi « 3,i -f |t| - t 


k » Ai 4- It ! t 


It also follows that Equation (2-78) can be expressed as 




ti*« i«0 


Since 
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2. 2. 1. 1. 2 Expansion of the Product (1 + px)° [1 + (y6/x)]*^ 

A straightforward application of the Binomial Theorem yields 


(;;)(:) /""x' 


k»0 msO . 


(2-87) 


The quantities t , p , and I are defined as above, with the exception that the 
ranges In this case are -n^t<n, |t|<p<2n, and 0£ I < (n- (It 1/2)] (where 
[ ] denotes the Integer part). Equation (2-87) can then be expressed as 
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The definition of the binomial coefficient 
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Also, 
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(n-ltl-i)l * (n-lV” (n- ltl)(n- Ul-i) (n- Itl- i+i) 


(n-UO! 


(-1) (itl-n) (itl vl) ••• (iti-n-lvi) 


which can be expressed using Pochhamer's symbol as 
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i-iiUii-n); 
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Therefoi'e, 
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(2-92) 


Similarly, 
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(2-93) 
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Finally, 
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Substituting this 


result into Equation (2-88) yields the result 


A/ V’ /n\ »tl \ ^-n 
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which can be 


expressed In terns ot the hypergeometrlo series ns 

i) * Z 1^) 


(2-97) 


•n.e expansions for (r/o.)*” in the true anomaly are obtained by substltmln* 
Equations (2-80) and (2-97) into Equations f2-75) and (2-76), respective y. w 


vields 
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ssing these expansion, explloltly in the true anomaly yields , 
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Expressing these 


expansions explicitly In the true anomaly yields 
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where 
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(2- 10 lb) 


2. 2. 1. 2 Expansion in the Eccentric Anomal\' 

It is apparent that if (r/a)*"and are expandible in a Fourier series in the 
eccentric anomaly, u , i. e. , 


i)'-I •: 


(2-102) 


Z»- 


(2-103) 


Although this expression is not of closed form, it is easily transformed to a 
closed-form expression using the linear transformation (Reference 26) 

F(a,b,c,x^) • (1- X*) F( c-fli, c-b, c, x*) 

which yields the result 

totP(p F(l-n, jS^) 
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Hence, the coefficients Wj * can be determined from the coefficients of the 
simpler expansions in Equations C-*102) and (2-103). 
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Another expression for the coefficients ‘ is provided by multiplying Equa- 
tion (2-103) by C and integrating over -.t<u<-t, which yields 
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In addition, the symmetry condition 
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follows, as before, by substituting the complex conjugates in Equation (2-107) 
anri comparing the results with those given in Equation (2-107). 

An explicit representation of the coefficients is developed next. This represen- 
tation is obtained directly by expanding the left-hand side of Equation (2-103) 
(rather than by using the expansions given by Equations (2-102) and (2-103) and 
then constructing the coefficients throug’. Equation (2-106)). The expansions in 
Equations (2-102) and (2-103) are also obtained, since they are special cases of 
Equation (2-105). 

Again, following Hansen's method, the definitions are made 
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it is easily verified that 

j) (2-111) 



Ll 

u 

u 

0 

I 


X a 3 




2-31 


(2-112) 


1 




u 


0 

u 


and 


P “ 


|t-s|-Ct-s) 

X 






(2-ll(5b) 


[] 


then U follows from Equation (2-93) that 
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Substituting Equation (2-119) Into Equation (2-115) and expressing the result In 
terms of the hypergeometrlc series yields 
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2. 2.1. 2. 2 Expansion of the Product (l-|9y) ^(l~/3y”^)^ 

Since n > I s 1 , It follows that 
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which can be expressed as 
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where oc and p retain the definitions given In Equations (2-116). Also, 
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and, similarly 
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Substituting Equations (2-125) and (2-126) Into (2-124) yields 
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Again, It follows from the definitions of a and ^ that 
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Thus, Equation (2-127) admits the hypergeometrlc series representation 
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Substituting Equations (2-122) and (2-129) Into Equations (2-113) and (2-114), 
respectively, yields the desired expansions 
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A closed-form expression can be obtained by a linear transformation (see 
Equation (2-10 lb) and the accompanying footnote) which >lelds the expression 
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2. 2. 1. 3 Expansion In the Mean Anomaly 

The Fourier expansion in the mean anomaly takes the form 



(2-132) 


It is clear from the discussion in Sections 2.2. 1. 1 and 2.2. 1.2 that the coefficients 
n s 

’ , referred to as Hansen's coefficients, can be expressed as 



(2-133) 


and also that the symmetry relation 



(2-134) 


holds. 

The above expansion has played a central role in the development of the classical 
planetary theories because of the desire for explicit time-dependent theories and 
because of the simple relationship between the mean anomaly and the time. 
Accordingly, this expansion and similar ones have been studied by a host of 
Investigators.^ 

There are several approaches which can be taken to develop an explicit repre- 
sentation for the Hansen coefficients. For example, the Hansen coefficients can 
be e.\pressed in terms of either of the previously derived coefficients, V” or 

^Leverrier (Reference 31) and Cayley (Reference 32) developed extensive tables 
for the expansions through the seventh power of the eccentricity. 
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The V* representation Is obtained from the Fourier expansion (Equa- 
t ^ 

tlon (2-65)) 
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(2-135) 


and a variation of the equation of the center (Reference 1) 
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Substituting Equation (2-136) Into Equation (2-135) and rearranging the summa- 
tion yields 




(2-137) 


% 

Comparison of this result with Equations (2-132) yields the relation 
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Similarly, the Hansen coefficients can also be expressed In terms of the coeffl- 
s 

dent ’ , The resulting expression of the Hansen coefficients has sometimes 
been x'eferred to as Hill’s formulation of the Hansen coefficients (Reference 33). 
However, this expression was given by Hansen (Reference 27) some 20 years 
before Hill. 

In addition to Hill's formulation of the Hansen coefficients, some of the repre- 
sentations obtained by Hansen and by Newcomb and Poincare will be presented. 
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2.2. 1.3. 1 Hill’s Representation for X 


n, s 


Hill (Reference 33) developed a representation of the form 
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where the function X * Is, to within a factor, a h^'pergeometric series, and 

t, p 


where 
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is the Bessel function of the first kind (Reference 30). 

This same form can be easl^"- constructed using Equations (2-130) and (2-131), 


l.e. , 
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and the expaasion 
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This expansion Is of major Importance In the expansion of the classical disturb- 
ing function (Reference 1). 
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where 


for r = 0 


Bq * I- e/a. for |r| =1 


0 for lr| > 1 


(2-143) 




and where 


(t ^ 0) 


(2-144) 


Substituting Equation (2-142) Into Equation (2-141) yields 
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and a comparison of Equations 02-132) and (2-145) gives the relation 


tn> tn.s r 

^ - L'''- 


(2-146) 


The range of r Is -n<r^n for n>0, -ooSr^oo for n<0 and -lSr<l for 
t = 0. 

The coefficients were shown In Section 2.2, 1.2 to be proportional to the 

Ixvpergeometrlc series; consequently, It Is clear that Equations (2-139) and (2-146) 
are of the same general form. 
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2.2. 1.3.2 Hansen’s Representation for X 


,.n,s 

t 


Another approach to the explicit development of Hansen's coefficients Is to expand 
the Integrand In Equation (2-133), Clearly, only the constant term la the expan- 
sion will remain after the definite Integral Is evaluated. Defining 




X a e 
'A 


ju. 
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(2- 147 a) 
(2- 147b) 
(2-147C) 


and substituting the relation 


dt -r-' 

Jl 


Into Equation (2-133) yields the contour Integral 


/ (l)" 


(2-148) 
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X £ di 
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where the contour c Is the unit circle | z | = 1 , 

Expansion of the Integrand In Equation (2-149) In powers of z will yield the results 
presented earlier In this section. However, the definite Integral formulation Is 
quite flexible In that the Integration variable can be transformed to either x or y 
via the relations 
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which follow from the well-known relations 
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Substituting Equation (2-150) into Equation (2-149) yields the contour integral 
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where the contour c is defined bv ( x | = 1 . 
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where 


/4 a k C0& <p 


Groupiog factors in the Integrand, developing expansions for each group, and then 
muWplylng all constituent expansions together yields the final expansion for the 
Integrand In Equation (2-WS), The form of the final expansion depends on hcnv 
the factors are grouped. Clearly, the product 








will yield a hypergeometric series representation as shown in 
Han?"n considers the factors 


Section 2.2, 1,2. 






The expansion for the exponential function vield 
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and, therefore. 






( 2 - 163 ) 
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It follows from the Binomial Theorem that 


i«0 


(2-164) 


Substituting Equation (2-164) into Equatioi (2-163) yields the result 
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If the following definition is made. 


p * m*l 


then the left-hand side of Equation (2-165) becomes 
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Thus, 


go 

^expj^] • 


C2-167a) 


where 


^ / .■ ^ (n*k4-mvl)! (p-m): «n! 


(2- 167b) 




■'[iii'Ti i I 'l. vi 1 i ; ^ 

0 
0 


f ^ , 

; i V ' 


U 

n 


uj 


L: 




L 


A 


The expansion of the product 


I i^-pA \ 


Is obtained by substituting -k for k, -jU for u, and x~^ for x In Equation (2-166). 
The result is 
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The product of Equations (2- 167a) and (2- 16Tb) yields the result 


(u^x) (ua«-) (1717-1;^) ; 
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since 


p = 




r-t 


it follows that r±t must be even and lt| . Hence, r Is defined as follows; 
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the right-hand side of Equation (2-169) takes the form 
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Multiplying this result by x yields the expansion of the integrand In Equa- 
tion (2.-159). 'Equation (2-159) then takes the form 
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(2-170) 
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Evaluating this Integral reduces to evaluating the integral 
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Clearly, 



( iTrj for t+s-k-1 = -1 
I 0 for t+s-k-lj^-1 


(2-173) 


from the theory of integration in complex variables or, equivalently, from Cauchy's 
Residue Theorem (Reference 30). 

In view of Equation (2-173), the expression for the Hansen coefficient in Equation 
(2-171) reduces to 
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where the coefficients AI and N are defined by Equations (2-167b) and (2-168b), 
respectively, 

Hansen obtained another expression for the Hansen coefficients by substituting 
Equation (2-151) into Equation (2-149) and using the relations 
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and Equation (2-156) to obtain the expression 
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Using the decomposition 
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where 
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Hansen develops expansions for the factors 


and 
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Substituting the expansions for these factors into Equation (2-177) and evaluating 
the result yields the expression 
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(2- 18 la) 


(2- 18 lb) 




-4 — — — 



It is Interesting to note that 


Gq = e 


-V 


Hq « e 


(2-182a) 
(2- 182b) 


Hansen provided other representation for the Hansen coefficients which can be 
found in Reference 27. 
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2. 2. 1.3,3 The Newcomb- Poincare Formulation of X * 

L 

Newcomb (Reierence 34) applied an operator approach to the problem of the ex- 
pansion of the classical disturbing function. This operator development relies 
on certain differential operators to produce an expansion In the eccentricity. 

The resulting development Is analagous to that obtained by using analytical ex- 
pressions for the Hansen coefficients and, consequently, provides another repre- 
sentation for the Hansen coefficients. 

In essence, the Newcomb operator method produces a power series In the square 
of the eccentricity, where the Newcomb operators are the coefficients. Evalua- 
tion of the Newcomb operators yields pure rational n.unbers. In addition, these 
coefficients can be evaluated recursively using the rec> 'uce relations that exist 
for the Newcomb operators. 

A complete discussion of the Newcomb operators would require an in-depth dis- 
cussion 'f the operator approach to the expansion of the classical disturbing 
function; this discussion Is beyond the scope of this section. However, In addi- 
tion to Newcomb's original work, the method is discussed and simplified by 
Poincare (Reference 35). Other treatments of the subject can be found in Refer- 
ences 1, 2, 29, and 36. 
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For the purpose of describing the Hansen coefficients in terms of the Newcomb 
operators, the following result given by Poincare (Reference 35) is considered; 


Os'oo msO 


(2-183) 


where the Newcomb operator 17^°^ ^^^(n|s) is a polynomial in n and s. Compar- 
ison of Poincare's result to a variation of the series of Hansen in Equation (2-132), 


= Ex- 


(2-184) 


vields the relation 
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(2-185) 


thus relating the Hansen coefficients to the Newcomb operators. 


The remainder of this discussion follows closely that of Iszak, et al. (Refeience 36). 
Inspection of Equation (2-183) Indicates that It can be expressed as 


iJ (^) * Z_. V 
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and X and z are defined b 5 ' Equations (2-147), 
It follows from the definitions of p and a that 
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According to Iszak, the change in the Indexes fi’om q, m to p,a simplifies the 
development of \'on Zelpel's recurrence relations (Reference 87) for the Newcomb 
operators. These recurrence relations follow from the partial differential equa- 
tion of \‘on Zelpel which Is derived In Appendix A and Is given by 
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Eliminating the explicit appearance of x in Equation (2-192) through the substi- 
tution 

* (t) 

and developing the resulting equation into a power series in e yields the partial 
differential equation 




n,krl 


+ (k-n)e^ii’x"’'‘''+ e* ae ^ h- 3i ^ 


X (2-195) 
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Substituting Equation (2-186) into Equation (2-195) and comparing similar terms 
vlelds the recurrence relation 
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The subscripts in this relation are restricted to nonnegative integers. 
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Another recurrence relation can be obtained from Equation (2-196) by interchang- 
ing the subscripts, changing k to -k, and using the symmetry relation 


^yO,<r 


(2-197) 


which follows from the symmetry relation for the Hansen coefficients (Equation 
(2-134)). The result is 


-4<rX. 
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Finally, a third recurrence relation can be obtained by summing Equations (2-196) 
and (2-198) to yield 

. n.k . , n.k+i ^ n,k-l 

4^y0+o*) X ^ 
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+ (k-n)X^ . ^ - (k+n)Xy, 
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thus eliminating the summation over r . 
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Initialization for the recurrence relation is provided by 

n.k 


( 2 - 200 ) 


and the fact that quantities with negative subscripts are treated as identically zero. 
Consequently, 


".1^ , VI 


(2-201) 


, n 


(2-202) 


etc. , are easily obtained. 

Although the Newcomb operators are rational numbers, the problem of generating 
them can be reduced to integer arithmetic by using the polynomials (Reference 36) 
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and the corresponding recurrence relations 
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and where 
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2, 2. 1.3,4 The Hansen Coefficient - A Special Case 

The Hansen coefficient ^ Is of major Importance In the development of the 
averaged equations of motion (In the absence of resonance phenomena), which 
are presented In Sections 3 and 4 of this document. Because of this importance 
and because It possesses a characteristic distinct from all other Hansen coeffl- 

n s 

dents (Iv ^ 0), It Is singled out for a special discussion. 


This particular Hansen coefficient Is the constant term in the Fourier expansion, 
In the mean anomal.v, of the product 


(i)' 


It possesses a finite h.vpergeometrlc series representation in either of the argu- 
ments or e^ , contrary to the general Hansen coefficients which have Infinite 
series representations. These finite representations can be obtained through a 
brute-force expansion of the Integrand of the expression 



(2-208) 


the special case of Equation (2-133) where k = 0 , However, this development 
la unnecessary since almost all of the work has already been performed in 
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Sections 2.2. 1, 1 and 2. 2. 1.2, More specifically, since 
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It then follows that changing the Integration variable In Equation (2-208) yields 
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( 2 - 210 ) 


which, In view of Equations (2-66) and (2-107), yields the relations 
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A Finite Hypergeometrlc Series Representation In 

It has already been shown that the coefficients Vg admit finite repre- 
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sentations In /3 and, therefore, they yield finite representation of the Hansen 


b 

U 

0 

0 

Li 

0 

0 

n 

U 

n 


a 

Q 

1 

I 

8 


2-56 


1 I 

iJ 



0 

L’ 

0 

li 

li 

U 


\ I 
^,5 

LI 

U 

0 

U 


y 

B 

I 

1 



t 









I 




V 





i 


coefficients through Equations (2-211) and (2-212). Replacing n by -(n-^l) 

in Equation (2-211) yields 


. -Cn-i) 

-(nrl\s Vs 

* T* 

0 coft 9 


(2-213) 


Substituting n-1 for n in Equation (2-lOlb) and dividing by cos 0 yields the ex- 
pression 


-(nti).s 

Xq * 


coi 


^ (T*^) 1 (2-214) 


where | s | in- 1. 

The expression for the coefficient X ^' ® is obtained by replacing n by n + 1 in 
Equation (2- 130b) and restricting the value of the subscript t to be identically 
zero. The result obtained is 
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This expression can be simplified in view of the fact that interchanging the first 
hvo arguments in the hypergeometric series has no effect on Its value. Since the 
only effect of a change in the sign of s on the hypergeometric series in Equa- 
tion (2-215) is a permutation of the first two arguments, it follows that 

F(li^ - r\-l, lsUl‘,y3^) « F(isl-n-l,-n*l (2-21G) 
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similarly, a change In the sign of s also causes a permutation in the binomial 
coefficients in Equation (2-215), Clearly, 


n-s 1 \ / 


n+\sl+l 


(2-217) 


In view of Equations (2-216) and (2-217), Equation (2-215) simplifies to^ 


(lt/3^) /3^) (2-218) 


A Finite H«’pergeometrtc Series Representation in e~ 

O 

A hypergeometrlc series representation in e“ can be obtained through a quad- 
ratic transformation of the hypergeometrlc series in Equations (2-214) and 
(2-218). Inspection of the transformed hypergeometrlc series in e“ indicates 
that they terminate after a finite number of terms. 

It is well known that any hypergeometrlc series F(a, b, c ; x) admits a quadra- 
tic transformation if and only if the quantities 

i(a-b), t(a^b-c) 


This expression is more easily obtained using the relation in Equation (2-211) 
and the closed-form expression for the coefficients (see Equation (2- 10 lb) 
and the accompanying footnote). However, it is also of value to pursue the 
simplifications which arise for the coefficients . 
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are such that any two of them are equal or one of tnem Is equal to 1/2 (Refer- 
ence 26). Inspection of the hypergeometrlc series in Equations (2*214) and 
(2-218) show that they satisfy the condition 


4-b a -(l-t) 


(2-219’ 


Consequently, they are Df the form 


F(a, b, jb^) 


(2-220) 


The quadratic transformation which yields the e* representation is^ (Refer- 
ence 26) 
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^There are two other quadratic transformations which can be applied in this 
case. They yield representations in terms of the arguments 
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with the simplification 
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which is easily verified from the definition of 0 . Applying Equation (2-221) to 
Equations (2-214) and (2-218) yields the expressions 
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An Associated legendre Polynomial Representation 

Anv hypei geometric series which admits a quadratic transformation can be ex- 
pressed as an associated Legendre Polynomial of the first kind (Reference 38). 
The hypergeometric series in Equations f2-214) and (2-218), which are of the 


form 


F(fl., b, a-ij+1; t) 
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These relations are particularly useful in view of the many recurrence relations 
available for the associated Legendre polynomials. The above expressions were 
also obtained by Cook (Reference 12), with the exception that a discrepancy of 
(-1)® appears between Equation (2-230) and the corresponding result of Cook. 

A review of Cook's work and the results he cites from Whitaker and Watson 
(Reference 30) Indicate that the missing factor does, in fact, belong in Cook's 
expression if Hobson's definition of the associated Legendre polynomials is used. 

Evaluation of Equations (2-229) and (2-230) for a few values of n and s yields 
the following Hansen coefficients : 
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‘i.2. 1. 4 Recurrence Relations for the Fourier Series Coefficients 

Recurrence relations prove very useful for the efficient evaluation of the coeffi- 
cients of the Fourier series expansions developed in the previous sections, 
provided the recurrence relations are sufficiently stable.^ Hansen gives several 
recurrence relations for the coefficients of the Fourier series expansions in the 
true, eccentric, and mean anomalies. More recently, Vlnh (Reference 39) has 
discussed recurrence relations for the coefficients ’ and ’ . Cook 

(Reference 12) and Cefola (Reference 11) have discussed recurrence relations 
for functions related to the special case of the Hansen coefficients Xq ’ . 

Glacaglla (References 15 and 40) also developed recurrence formulas for the 

Hansen coefficients, and Cefola (Reference 41) has developed a new recurrence 

n s 

relation for the coefficients * . 

iD g g 

Previousl}', the coefficients ’ , and X^ ’ were shown to be expressed 

simply in terms of the hypergeometric scries. It follows that Gauss' contiguous 
relations (Reference 26) for the lu'pergeometric series could serve as the basis 
for constructing recurrence relations for these functions. However, recurrence 
relations that pei’mit only one varying index are often the most desirable. Since 
the parameters of the hypergeometric series (a, b, c) (Equation (2-28) are linear 
combinations of the indexes n, s, and sometimes t, the form of Gauss' contig- 
uous relations is not well suited for developing recurrence relations with a single 
varying index. 

Orthogonal polynomials are suited for single-varying-parameter recurrence rela- 
tions. Furthermore, it was shown in the previous section that the coefficients 
s 

Xq * are simply related to the associated Legendre polynomials. Consequently, 

^Slnce truncation and x'oundlng errors are Introduced into the evaluation of the 
recurrence relations, it is important to know how these errors are propagated 
through the recurrence process. II the errors do not grow relative to the mag- 
nitude of the function being evaluated, the recurrence relation is said to be 
stable; otherwise, .the recurrence relation is unstable. 
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the recurrence relations for this set of orthogi al polynomials provide the foun- 

-n s 

elation for the recurrence relations for the coefficients X. ’ and for the coeffi- 
dents V~ and , in view of Equations (2-211) and (2-212). 

-n s S 

Recurrence relations for the coefficients ’ and ’ are obtained through 
a more classical approach used by Hansen, 

2.2.1. 4,1 Recurrence Relations for the Coefficients X^*®, V^, 

Recurrence relations for these coefficients are obtained from the following re- 
currence relations for the associated Legendre polynomials obtained from Refer- 
ences 26 and 42. The fixed-order recurrence relation is given by 
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the fLxed-degree recurrence relation is 
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and the varying-order-and-degree relation is 
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Equations (2-232) and (2-235) can be combined to yield the recurrence relation 
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It is easlh- demonstrated from the results in Reference 26 that 
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and It follows fi'om the principle of induction that 
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In addition, It Is easily shown from the definition of the associated Legendre 
polynomial that 
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Combining Equations (2-237) and (2-238) yields 
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and it follows from Equations (2-238) and (2-239) that 


(2-239) 


. un.q p; w 


(2-240) 


Inverting Equation C--229) and substituting the result into Equations (2-232) 
through (2-233) yields the following recurrence relations: 


L 
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) 


-(intO,s+l t r 9,s 


nt s 




(2-242) 


^0 = U.s-iUmt-Z ) J (2-243) 


(ni-OjS 


^ [(«-i) x;"'% (n-») e X'/”*''’ ] (2-244, 


(a + S-I) 


In the above equations, the superscript s is restricted to nonnegative values, 
i. e. , s > 0 , This restriction is quite satisfactory in view of the symmetry re- 
lation (Equation (2-134)) 


X 


rt.S 

k 


X 


A.-S 

-k 


In view of Equations (2-229) and (2-236), it follows that 


-(AtI), n-1 
aq 



and it follows from Equations (2-229) and (2-239) that 


X 


0 


I e -(A+i\n-i 


Also, Equations (2-237) and (2-229) can be combined to yield 




(2-245) 


(2-246) 


(2-247) 
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In addition, inspection of the definition given by Equation (2-229) indicates that 


s 0 


(2-248) 


and, more generally. 


2 0 


(2-249) 


for all s such that ls| > n . 

Tlie recurrence relations for ® are obtained by substituting Equations (2-230) 
into Equations (2-232), which yields 


n.s ^ 


(an+OXo - 




n-3^,. 

0 


(2-250) 


I r as 


(2-251) 


n.s 1 a^) n-l,S 

0 * nMi •! ^0 


0 


(ivii-OeXQ ' (2-252) 


^o‘* ex; 


n.5-1 


(2-253) 


The superscript s Is again restrict 'd to nonnegative values. 
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Applying Equations (2-236) and (2-238) to Equation (2-230) yields the result for 
the special case 


ol 


1 1 


- I i 


! ! 


' rl 

I 


^ (n+lV. 


(2-254) 


Also, Inverting Equation (2-230) and substituting the result into Equation (2-240) 
vields the recurrence relation 


,n.n n-i 

■0 " ‘^0 


(2-255) 


and two successive applications of Equation (2-240) with Equation (2-230) yield 


yi.kl ( oLn t iD ( fllyi” 0 n. n*5,rt'3L 

- n 


(2-256) 


Clearly, 


n,wti 

X, * 0 


follows from Equation (2-230) and, generally. 


X(J a 0 


(2-257a) 


(2-257b) 


for s such that | s | > n+1 . 

Inspection of the recurrence I'elatlons in Equations (2-242) and (2-251) indicate 
the appearance of e as a divisor. Consequently, these recurrence relations 
appear to be of little value for cases with small eccentricities. This difficulty 
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Is easily avoided by using the expressions 

S A s 


« e*A 


n 


s s 

where the definitions of the functions and follow from Equations (2-224) 

and (2-223), respectively. The resulting recurrence relations for the functions 
s s 

and B , are free of the eccentricity divisor, 
n n+1 

— Q S 

The recurrence relations for the coefficients V” and W“ * are easily ob- 

S t) 

talned by substituting Equations (2-211), (2-212), and (2-213) Into the recurrence 
relations for the special case of the Hansen coefficients and ^ 

2.2. 1.4.2 Recurrence Relations for the General Hansen Coefficients 

n s 

since the general Hansen coefficients Xj^ ’ do not apparently admit a simple 
orthogonal polynomial representation, the classical approach of Hansen will be 
used to develop the recurrence relations for these coefficients. 


For convenience, the following definitions are made; 

r 

* I 


X » fc 


0. 


Then, 


z » e 


. rt S n-l i ^ 


n s-t 


(2-258) 


Using the relations 


I t| M'i ■ -j , ., 

J I / L. -i _ -1- L. -I— k — i — L 


-*i>C ^.4 . 




-. i 4 

d\ i «'^ 

dp ^ /^ \ 

dx. i ol^ cos*^ ' 


(2-25 9b) 


(2-259C) 


Equation (2-237) can be simplified to yield 


n s n ^n-l , .s-l 

flii ^ 2COS0 


n-x / s-i ■»+x\ I n*l 4 

i T- .-0 X' a T . -^ T p (X - X ) + S COS9yO X 


Since ^ 




cos ^ 




s 0 


or, equivalently. 




Sim ^ (xrx"') - leoi ^ X - 0 


then the products of Equation (2-262) with the factors 


n-JL 5 n-i % 

^ Y\p K S/) X 

t — *- — ■; — and — ' 

acos<^ 2ws9 


In classical elements, this equation takes on the more familiar form 


^ 1 + e C04 ? 


a 0 


C2-260) 


(2-261) 


(2-262) 
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are identically zero and can therefore be added to the right-hand side of Equa- 
tion (2-260) to yield 


. n-1 s-l , V j n-A s n n*l $ 

2 “ /3 A ■ ’’ *~T +* (s-n^ CosO/5 ;< » ■ o X (2-263) 
dl ' COS0 ' ' ' cos^ ~ 


d 0 i 
* * 


n n-l srl , ^ ^ o-a s -n n-i & 


1:U 


d rt a (or s") n-1 s-1 (o-s^ n-i «-»l S n-i S 


*51'"* * 


Also, since 




am ^ /-IN yO 

O (.X 4- X ) ■♦• ”£ —• a i 

a COS 0 COS 


(2-266) 


the product of Equation (2-266) with the right-hand side of Equation (2-265) .yields 


da/ 4 cos®^ / Acoa^<^ ' 


as ^ s sin n s (n- as ) sinct n s*i 


2 

( n- s ') iin 


(2-267) 
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A less familiar form of the expression 1 r 6 Cost * 

ad- «•*) 
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Hansea also obtains the result 


i .n 




- n(3in-‘3)p” ^ X* <► n(w-l')|0*' 
*► s(n-i') 


(2-268) 


by differentiating Equations (2-263) and (2-264) with the operator € 7 — and sub- 

dl 

stltutlng those equations back Into the result. Equations (2-262) are also helpful 
In obtaining the final form. In addition, multiplying Equation (2-262) by the 
factor 

s(n-l)^ X 

and alternately adding and subtracting the result from Equation (2-268) yields 
the expressions 


a. s d n s 

^ T-T fi ^ +1— ,0 X = 
di^ ' dt • 


- ^ n (3n*i) ~ 1 s (v\»l) j yO x 


(2-269) 


-► n(n- ^x*4- lsCn-0 sin^ yo” ^x* ^ 


1 d^ n s d ns 


» [n(n-a')+ 5^i*as(n-l)] cos^4 yO*' ** x^ 
• [n(ln*3)<t* 3L%(n*l)] /j" 


(2-270) 


•1 {a-1) X* - iiCn-O ''x 
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n-4,ft 


^tr\(n-a)+ s^ j X^ ’ - v\(an-5)X^ ' + n(n-l)X^ ' (2-277) 

+ 5 ( 11 - 1 .) ti»T^ X^ ^ - s(n-l) * 0 


j^w('l- i.) + ^^'-is(ii-l) jcos^^X^ ‘ - j^n(an-'5)« Sls(irt-l') jx^ 


(2-278) 


[n (w-l) + 5^+- 3.&(n-l)j cos^^ X^ * - [n(in-a) + 3.s(n-l) jx^ 
•*• n(w-l) - 3l&(n-l) sin^ X^ ’ - t^X^* = 0 


(2-279) 


Setting Equations (2-260) and (2-263) equal and substituting Equation (2-27 la) into 
the result yields the recurrence relation 


n-l 4-1 n-i.4*l n»l,s a , 

Sin 0 X^ 4* si») 0X^ '♦•«LX^ - ^ C04 0 X^ (2—280) 
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Also, replacing n with n-2 in Equation (2-272) and setting the result equal to 
Equation (2-277) yields the recurrence relation 


n-3.» 


n [(n-a)'- s^]c6S^0X^ ‘ - n(n-a)(an-3)X^ ^ 


(n-l) [n(n-a) 4* atsCOS0 ] x!J t*‘(n-a) x'^*^ * 0 


(2-281) 


This recurrence relation is particularly attractive since the superscript s is 
fixed and since the eccentricity e = sin0 does not appear. 
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Several more recurrence I'elatlons can be obtained by combining the above recur- 
rence relations or, equivalently, the differential equations. All of the above 
recurrence relations are valid for positive and negative values of n. 

2.2.1. 4, 3 Recurrence Relations for the Coefficients 

It 

The procedure for generating the recurrence relations for this case is Identical 
to the procedure used above. Since 




(2-282) 


.• U' 


dy U-;dy')^ 


(2-283) 


where 


it follows that 




dp n 


= n,3e«^4/5"'‘ - l) * (2-284) 




which can be expressed as 


[y(U)3^)- - »y) + a (2-285) 

ay 


Substituting the series representations 




(2-286a) 


2-75 


1 ---^1 - 1 : 1 , 4 ..:. .-i .1 ■ ' 


dp 
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yields the recurrence relation 
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n.s VI. 


n*-l s.i 
^ X ^ ct>s 


(2-286b) 


(t+ 1 • ^t-S -»• Ct+*s)j3^ jW^ ’ V * 0 (2-287) 


Another recurrence relation can be obtained from the expression 


«-288) 
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which Is easily verified. Substituting Equation (2-286a) Into Equation (2-288) 
yields the recurrence relation 

nvl,s-t X d r n.& a n,s i 


u 


In addition, Hansen also gives the recurrence relation 


4 r. I*'-® t”'® 1 


(2-290) 


Recently, Cefola (Reference 41), using the technique of Hansen, obtained the re- 
currence relation In which the superscript s Is fixed, l.e., 


(a«-l) tscos^ - n(n-0]W^ 


n-l,s 


(2-291) 
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Thts recuri'ence relation Is a generalisation of the recurrence formula given by 
Vlnh In Equation (3.8) of Ueference 39. 

Cefola's recurrence relation Is reminiscent of recurrence relations for orthog- 
onal polynomials and thus leads to the conjecture that the coefficients w”* ^ may 
have an orthogonal (x)lynomlal representation. This conjecture Is shown to be 
true In Appendix B of this document. 

2.2.2 Fourier Expansions In the Longitudes 

Tl\e results obtained In Section 2. 2. 1 can now be used to construct the Fourier 
series expansions In the true, eccentric, and mean longitudes of the product 
(Equation (2-54)), 





(2-292) 


where X is ai\v of tlie three longltiules. it was previously shown that the coeffl- 
11 s 

dents Aj’" of the expansion in a particular longltiule are related to the coeffl- 
tl s 

dents a^ ’ of the expansion In the corresponding anomaly by the expression 


(k-jnh) 


Is-tl 


(2-293) 


where >7 = slgn(t-s). Equation (2-293) Is equivalent to Equations (2-(52a) and 
C2-d2b), The form of Equation C-~293) Is used to avoid reciprocals of the complex 
polynomial (k^-jh)”' by replacing It with Its conjugate iwlynomlal, 1. e. , 





(2-294) 


It appears that the form of Equation (2-293) admits a "computational singularity" 

n s 

for vanishing eccentricity; however, the coefficients a * contain a factor that 

(i 

absorbs the eccentricity divisor, as will be demonstrated. 
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2.2,2. 1 Expansion In the True Longitude 

Since only the expansion of the parallax factor, r/a , Is required, It Is considered 
first. In view of Equations (2-100), (2-292), and (2-293), It follows that 




-Ul. . , ,tn jiL 


(2-295) 


whe;e -Qo<t<oo for +n and -n<t<n for -n, and where 7i~ slgnt, 


In 


Inspection of the expressions for the coefficients (Equations (2- 100b) and 
(2- 10 lb)) show that they are proportional to /3*^‘ or, equivalently, e**"^ (Equa- 
tion (2-73)), Hence, If the coefficient p Is defined as 


in It) tn 

* e 


(2-296) 


then Equations (2-295) take the form 




(2-297) 


IsL 

Multiplying both sides of Equation (2-295) by yields the final result 

I r\~^ jiL . .l-tl xn jCitS^L 

[i) ^ " 

t 


(2-298) 


where -oo<t<oo for +n and -nStSn for -n. 

2. 2. 2. 2 Expansion In the Eccentric Longitude 

In view of Equations (2-130), (2-131), and (2-293), It follows that 


ia 


a' ‘ 


I 


e - Ck-^T^K) W^. e 


(2-299) 
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where r? = sign (t- s), -n<t<n for -^n and -oo<t<oo for -n. Inspection of 

+n s 

Equations (2-130b> and (2- 13 lb) show that the coefficients W^' * are proportional 

^ ^|t-s 1 It-s I , ,, in, 3 . , , , 

to /3 or e ; hence, If w^. is defined by 


tn,s It'Sl tn.s 


e w* 


(2-300) 


then Equation (2-299) takes the form 


/r\" jsL V., . I 


\t-tl tn,i 

w. e 


(2-301) 


where the definition of rj and the ranges of t are the same as In Equation (2-299). 
2. 2. 2. 3 Expansion In the Mean Ijongltude 


The expansion In the mean longitude follows from Equations (2-132) and (2-293) 


to vleld 


lr\-^ jsL j-' 

\l) ^ ^ e 


(2-302) 


Also, Inspection of Equation (2-174), (2-180), or (2-185) shows that the Hansen 
coefficient Is proportional to ; hence, If the coefficient 


Is defined bv 


it-sl tn.a 

a e 


(2-303) 


then Equation (2-242) takes the form 
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Z . vU-t\ t< 


U-t\ tn.a 

K. e 


where rj ~ sign (t - s) 


(2-304) 
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2. 2. 2. 4 Recui’rence Relations 

Recurrence relations for the coefficients of the expansions in the longitudes 
follow Immediate l.v from the recurrence relations for the coefficients of the ex- 
pansions in the anomalies presented In Section 2,2. 1,4. For example, substi- 
tuting Equation (2-303) into the recurrence relations for the general Hansen 
coefficients yields the recurrence relations for the functions K" * , Combin- 

Ing these recurrence relations with the simple recurrence properties of the 

I s-t I 

complex polynomial (k - jrjh) yields the recurrence relation of the product 


tn.s 

Y 







(2-305) 


As an example, substituting Equation (2-303) into Equation (2-281) yields the 
result 

[ a 3 1 A / V n-3,4 

(n-ar- j cos ^ Y. - n(n-a')(3»'-3W. 

■* (2-306) 

e (n-i) l^ia(n-a) 3lstcos(J> j yj* Y » 0 
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SECTION 3 - EXPLICIT THEORY FOR THE NONSPHERICAL 
GHAVITATIONAL PERTURBATION 



This section presents the explicit development of the first-order averaged equa- 
tions of motion for the nonspherical gravitational iJerturbation. Section 3, 1 
presents a general discussion of the development in spherical coordinates of the 
nonspherical gravitational disturbing function. In Section 3.2, the disturbing 
function is developed explicitly in terms of the equinoctial elements. Specifically, 
Section 3.2. 1 presents the rotation of the disturbing to the equinoctial frame and 
Section 3. 2. 2 Introduces the necessary Fourier series expansions. 

In Section 3.3, the averaged disturbing function is obtained. The concepts of 
time-dependent and time-independent averaging are introduced and compared. 

It is shown that time-dependent averaging does not always remove all short-period 
terms in the disturbing functions and that time-independent averaging may exag- 
gerate the amplitudes of the remaining medium- and long-period terms in the 
averaged disturbing function for the cases of nonresonant and resonant tesseral 
harmonic terms, respectively. In addition, the zonal harmonic, combined zonal 
and nonresonant tesseral harmonic, and resonant tesseral harmonic disturbing 
functions are Isolated. 

Section 3.4 develops the partial derivatives necessary for the a\’eraged equations 
of motion. The equations of motion ai*e developed separately for the zonal har- 
monic, the combined zonal and nonresonant tesseral harmonic, and the resonant 
tesseral harmonic models. 
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3. 1 THE NONSPHERICAL GRAVITATIONAL DISTURBING FUNCTION 

It Is well known that the gravitational force, can be represented as the gra- 
dient of a potential function \’, i,e. , 

F s -VV(x,4,a) (3.1) 


where the gradient operator V is defined in Cartesian coordinates as 


V 


6 a. 3 a 3 a 

1“ ' r" j 1“ ^ 

dx dy bi 


(3-2) 


AAA 

where (i, j, k) is the orthogonal triad of the Cai'tesian reference system. 


The particular form of the potential function associated with the gravitational 
force exerted by the attracting body depends on the mass distribution of that 
body. The potential function must satisfy Poisson's equation 




(3-3) 


where the divergence operator is defined by 


a ^ 




(3-4) 


and where G is the universal gravitational constant and ^(x, y, z) is the density 
per unit volume at the point (x, y, z) . At all points where the density vanishes, 
i.e., outside the attracting body, Poisson's equation reduces to Laplace's equa- 
tion, i.e. , 


V 




(3-5) 
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The general solution of Laplace's equation yields the potential function for the 
gravitational force exerted by a body of arbitrary mass configuration on an ex- 
terior particle located at the position (x, y, z). The potential function for a given 
mass configuration Is specified by the appropriate boundary conditions In addition 
to the general solution of Laplace's equation. In general, these boundary condi- 
tions are unknown, and. In practice, the potential function is ultimately deter- 
mined by a semlempirical method. This method assumes knowledge of the form 
of the general solution of Laplace's equation which is obtained below. 

The method of solution for Laplace's equations is the standard separation of 
variables technique and is usually developed in spherical coordinates (r, ,0 ) 

where r>0, 0<</'<2«', and -7r/2^0^ n’/2 . In spherical coordinates, Laplace's 


equation takes the form 


0. 

6r \ / "** COS 


1 

COS0 d0 


dv\ 1 

df »■ 


The solution Is assumed to be of the form 

V = AM vm $W) ( 

Substituting Equation (3-7) Into Equation (3-6) and dividing the result by Equa> 
tlon (3-7) yields the differential equation 


ACr) dr \ dr / $ cos0 d(J> \ 




d’vS') 




3-3 


since only the first term Is dependent on r , It must be (most generally) a con- 
stant to satisfy the above differential equation. It Is convenient to choose the 
constant to be of the form 1) , hence 


— — 

Air) dr \ 


dAir) 

dr 


) = Kui) 


(3-y 


which can be expressed In the form 


j, d^A(r) 

'* dr^ 


Ir 


dA(r) 

dr 


-l(Ui)A(r) = 0 


(3-10) 


The agreement between the power of r and the order of the derivative in each 
term of this equation suggests a solution of the form 


A(r) a 


(3-11) 


Substituting this form Into Equation (3-10) yields the equation 


n(n-l) •♦•In - = 0 


(3-12) 


which admits the solutions 


I ^ 


(3- 13a) 
(3- 13b) 


Therefore, the general solution of Equations (3-10) is of the form 


A(r) a ♦ Cg^p’'*"'^ 


(3-14) 
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where c- and are arbitrary constants. For the gravitational potential, 


s 0 


since the gravitational potential is assumed to vanish at infinity, i.e. , 


(3-15) 


Jim V(r,6,(^) s 0 

r J»«D 


(3-16) 


The remainder of the potential function is then determined. Substituting Equa- 
tion (3-9) into Equation (3-8) and multiplying the result by cos“<f> yields the dif- 
ferential equation 


•4 r’l 


^ i -^V • 0 


(3-17) 


The last term is clearly constant since it alone depends on t/t , If this constant 
is denoted by -m" , then 


a. 

^ • 0 


(3-18) 


which is the equation for a harmonic oscillator and which admits the solution 




(3-19) 


where C and S designate arbitrary constants of Integration 
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The final function ^{<p) Is detei*mlned by substituting the constant -m“ Into 
Equation (3-lTl and multiplying the I’esult by $/ cos-^ which yields the differ- 
ential equation 


1 d I \ 

X TT COS0 TT ♦ I Ul*i) jx 

cos 9 d<^ \ d<p I I cos^9 


$ a 0 


(3-20) 
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Transforming the Independent variable In Equation (3-20) through the relation 


X 3 


(3-21) 
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yields the resulting differential equation 
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l*x‘ 


-I 


(3-22) 
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For m = 0 , this eciuatlon reduces to the classical equation named after Legendre 
(Reference 30) and which admits as solutions the Legendre polynomials P^fx), l.e. , 


r t 

I I 


• CjP^Cx^ » (form = 0) (3-23) 
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where Cg Is an arbitrary constant. The Legendre polynomial Is defined on the 
Interval -l^xSO by Rodrigues’ formula (Reference 30) 




(3-24) 
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The solutloa to the more general differential equation given In Equation (3-22) Is 


denoted by 


^ • t5%CVw'J') 


(3-25) 


and Is called the associated Legendre polynomial of degree X and order m . The 
associated Legendre polynomial Is defined on the Interval -l<x< 1 by Ferrer 


(Reference 30) to be"^ 


'i.m ^ - Vi. X ; 


(3-26) 


A complete discussion of the general differential equation Is given In Reference 30. 
However, one fir ther point Is of Interest for this discussion. The general differen- 
tial equation ^iven In Equation (3-22) can be obtained from Legendre’s differential 
equation lor which m = 0 by differentiating the latter m times. However, since 
Legendre's differential equation admits a polynomial solution of degree X, m 
must be restricted to the range 0 < m^ X In order to obtain a nontrivial result. 

In view of Equation (3-14), (3-15), (3-19), and (3-25), the potential V depending 
on the constants X and m and denoted by can be e.xpressed as 




(3-27) 


Some authors Include the factor (-I)**' In the definition of the associated Legendre 
polynomial given by Equation (3-26) (Reference 42). Tliese differing definitions 
have contributed to a certain amount of unnecessary confusion. The effects of 
this definition difference will be observed In the sign of the spherical harmonic 
coefficients of odd order m. Reference 30 adopts the notation of for 

Ferrer's definition and the notation ’ 

- i-i) 

for the alternate definition. This convention Is observed In this report. 
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where 


Ci,™ - «.C3C 
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c^c S 


(3-28a) 

(3-28b) 


are dimensional spherical harmonic coefficients. Since the dimensions of a grav- 
itational potential function are (length/tlme)~ , a set of dimensionless spherical 

harmonic coefficients (C« S. ) are obtained through the definitions 

4t,in 


m 




(3-29b) 


where ju Is the gravitational parameter of the attracting body, defined In terms 
of the 1 
G, by^ 


of the mass of the attracting body, M , and the universal gravitational constant, 
1 


ytc « GM 


(3-30) 


and designates the mean equatorial radius of the central body. 

A complete account of the method Is given by FU/.patrlck (Reference 43) and a 
more physical approach Is provided by Battln (Reference 44). 

The complete general solution Is obtained by summing over all admissible values 
of X and m , which obtains the following final result: 


oo X 

V , - ^ y 7 V. 

fafl W«0 


(3-31) 
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\'alues for the gravitational parameter are better obtained through observation of 
the satellite mean motion, n, and the semimajor axis of the satellite orbit and 
through the use of Kepler's Third Law expressed as n“ a^ - yu. 
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and where ^ = the gravitational parameter = GM 

r = the distance of the satellite from the origin of the 
coordinate system reference frame 

= the mean equatorial radius of the attracting body 

0 = the latitude of the satellite 

0 = the body-fL\ed longitude of the satellite 

C- , S, = the spherical harmonic coefficients which are deter- 
mined empirically for a given body 

m^^ ~ associated Legendre polynomial of degree X , 

’ order m , and argument x 

The term of zero degree and zero order, l.e. , 


(3 -33) 


corrcs|>onds to the potential function of the point mass of classical two-body theory. 

Inspection of the expression In Equation (3-32) Indicates that the values for the 
spherical harmonic coefficients ^ are completel.v arbitrary and can be taken 
to be Identically zei'o, l.e., 


S, « a 0 
1,0 











Furthermore, It can be shown that the spherical harmonic coefficients C . , 

C , and S vanish Identically If the origin of the coordinate system Is placed 

If 1 if 1 

at the center of mass of the attracting body (lleforence 43). Consequently, under 
this condition, 

Vj 0 5 0 (3-35i 




(3-35a) 
(3 -35b) 


and the nonspherlcal gravitational potential function takes the form 
^ It 4. m«0 

The disturbing function R is obtained by taking tlir negative of the potential func- 
tion and deleting the point-mass term to yield 

00 1 ^ 

R * 7 ^ X! ^ ^1.01 ^ 

n\tO 

Tlie partial derivatives with respect to the cKiulnoctlal elements of the disturbing 
function are requli'ed by the equations of motion. It Is desirable to express the 
disturbing function directly In terms of the equinoctial elements, rather than 
relying on the application of the chain rule. A complex variable representation 
of the disturbing function will facilitate the transformation to equinoctial ele- 
ments. The disturbing function Is given by the real part of the expression^ 
ce X 


R* • T ^ ^ - j s,„„) e 


(3-38) 


i*X *w«0 


Reference to the real part of the complex variable representation will be dropiwd 
until the final form of the iwtentlal and Its partial derivatives are obtained. 
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1 



where j Is the Imaginary unit, v'-l . nnd the boih -fixed longitude, ^ , has been 
expressed as the difference of the right aseenslon of the satellite, g , and the 
Greenwich Hour Angle or some txtuivalent angle for a central body other than the 



3.2 TRANSFORAUTION TO THE EQUINOCTLAL ELEMENTS 

All quantities In the disturbing function which are dependent on the satellite pos- 
ition must be expressed in terms of the equinoctial elements. This transforma- 
tlo .nplicitly or explicitly requires a rotation of the coordinate reference frame 
associated with the coordinates to one of the equinoctial reference frames 

(direct or retrograde). A general discussion of rotations is presented in Sec- 
tion 2. 1. 1 and the specific rotation required is discussed in AppendLx A of Refer- 
ence S. 

Inspection of Equation (3-3S) indicates that the only quantities dependent on the 
satellite position are the spherical harmonic functions 






e 




Tlie general theory of the rotation of the spherical harmonic functions Is dis- 
cussed in Section 2, 1.2. 

In addition, the radial distance of the satellite, which is Invariant under a rota- 
tion, must be expressed in terms of the equinoctial elements. This is accom- 
plished through a Fourier series representation in the true, eccentric, or mean 
longitude. Since sines and cosines of multiples of the true longitude, L, are 
Introduced by the rotation of the coordinate reference system, it is really neces- 
sary to develop Fourier series expansions of functions of the form 

(x) «>• i (r) 

The development of these Fourier series expansions is discussed in Section 2.2. 
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a. 2.1 notation of the Snherieal Uarmonlo Funetlons 


The rotation to the equinoctial reference frame can be exprcaseil by the prixluct 
of three simple rotation matrices (see Apiiemllx A of Uoference 5) 




(3-3in 


'■ u 

i 0 

;i 

•: I u 


u 


where A ami l are the rijjht ascension of the ascemUng notle and the inclination 
of the satellite orbit, resix>ctiveVv, and I is the reti'ograde factor, *l1\e rotation 
matrices Rj and R^ are defined by Equations and (2-7). The rotation 
transforms the spherical coordinates (r, ) to the spherical coordinates 

(V, L, <>'i 0) , The satellite latitude is lilentically zero in the equinoctial refer- 
ence system since the fumlamental plane is the orbital plane. 

There are two approaches to the transformation of the spherical harmonic func- 
tions, Kaula (Reference 16) usetl a brute-force substitution of the Keplerian 
elements for the right ascension and latitude of the satellite.^ topressions were 
obtalneit of the form (in complex variables) 


! U 


i “ 

ill 




(3-40) 


where f is the true anomaly ami the Inclination function, r, (1) is defined to 

A.m.p' 


tU-aO! vn t 


r-’ tu-ao! vn 

F. ti^ • N 

tlU-05 


VH. 




(3-41) 


This substitvrtlon can be obtaliunl from Equations (2-22) ami (2-23) by ileslgnatlng 
A as the right ascension and X* as the true anomaly- , l.e. , A' =■ f. 







where k is the integer part of ( i - m)/2 , the summation index, t , varies through 
the range 0 < t < min(p,k) , and c is summed over all values for which the binom- 
ial coefficients are defined. 


Iszak (Reference 45) and Allan (Reference 46) simplified the expression for the 
inclination function by using the inclination half angle to obtain^ 




ii-W- ap-Sk 


(3-42) 


(\/a) 


m-X>ap 


where k^ = max(0, X -m-2p) 
kg = min(2X-2p, X-m) 

To facilitate the development of a simple recursive scheme for evaluating the 
inclination functions, Cefola (Reference 13) suggested the alternate approach 
based on the theory of the rotation of the spherical harmonic functions presented 
in Section 2. The rotation of the spherical harmonic functions to the equinoctial 
reference frame takes the form (Reference 25) 




X 

ZL. (l-m)! 

s«-X 


Pi,, 'O' 


(p,<T,T) t‘ 


isL 


(3-43) 


This particular expi'essicn was ubt^iiped by Iszak} Allan's expression incoi’por- 
ates the factor in the above definition. 
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since the latitude of the satellite In the LMUtuoettal reference frame, Is Iden- 
tlealVv r.ero. I'he parameters p, 0*. and T for this ease will be defined shortly. 

It should bo noted that the function P (0) Is a constant depending only on the 
Indexes .i and s. 


It follows from Equation (S. G. 11 of Reference 26 and from the footnote on page S-7 
that 






n 

1 1 

r’ 

r 1 


(3-44) 


(nils expression Is valid for all Integers s . ) 
it Is obvious that for odd values of X. * s , 

P , (0) 3 0 (3-45) 


I'slng the relations 



Un-Qi! 



(:i-46a) 




C3-46b) 


It Is easily shown that 


P. (0^ 


(-0 
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provided the definition 


(-1)!! St 


Is made. An alternate definition which avoids the double factorial notation is 


P , ,( 0 ) * (-0 


u-s)y3L 




(¥)'• (¥)> 


(3-48) 


Since only terms which satisfy the condition .t- s = 2p contribute to the summa- 
tion in Equation (3-43), the range -ASs^Jt can be replaced by the range 
0 <p< £ to j'leld 






(3-49) 




This modification is of no significance for machine processing and will not be 
adopted in this report. However, it does demonstrate the relationship between 
the present formulation and the standard Kaula approach given by Equation (3-40). 

3. 2. 1. 1 Determination of the Function ) 

in s 

To obtain the appx’oprlate fox'm of the function S.^^’ (^,<r,T) , it is necessary 
to determine the parameters p, CT, and r in terms of the Euler angles (A', i', o»), 
which describe the orientation of the equatorial reference system relative to the 
equinoctial reference systems. The rotation from the equatorial to the equinoc- 
tial reference systems, which was derived in AppendLv A of Reference 5, is given 
by Equation (3rl3). 


i 
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Tor the purpose of this discussion, the Inverse ti'ansformatlon Is requii'ed, l.e. , 




Since the constituent rotation matrices, Rj^(6) and R^(9)t are orthogonal. It 
follows that 


* R^(6) - R(-8) 


Consequently, 


T‘^ - RjC-Il) Rj(-i) Rjdm 


Comparison of this result wltli Equation (--5) indicates that, for this case, the 


Euler angles are 


0)' » -n 


a * la 


The primed symbols denote the general Euler angles and the unprlmed symbols 
denote the familiar Keplerlan orbital elements. The parameters p,CT, und T 
are thus obtained from Equations Oi-lt*)* l.e.. 


■ ■ m — — 1— — i* 
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(3-54b) 


(3-54c) 


These expressions can be substituted into Equations (2-44) to obtain the function 
* however, It Is just as easy to substitute the above Euler angles directly 
into Equation (2-45), which yields the result 

exp [jU-m)5 ] 


-m-s m-s 

^i/a. 


(-ASs^-m) (3-55a) 


S-,/ 3 , Pj.^ Ceo (-mis<m) • (3-55b) 


Cill S;,^ Pj„ CCO 


(n^sSJL) (3-ooc) 


The substitutions 


c. - c, 


S = -S 


were used to obtain the above results. 
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The collapsed form of the iaclidatlon function is obtained by substituting the Euler 
angles (Equations (3-54)) into Equation (2-50) to yield 


] €)tp ] & 


(Us^Ul-O! 


$«4rn 

(*l) ^\ix ^i/a ^-4 (^C-0 


(0 i s ^ m) (3-56a) 


(m<s<X) (3-56b) 


where e assumes the values fe - i 1 , 


The inclination functions can be expressed in terms of either the equinoctial 
elements p and q or the direction cosines, with respect to the equinoctial i*ef- 
erence system, of the i axis of the equatorial system. The relationship betix'ccr. 
the direction cosines («,j3. V). the elements p and q, and the Keplerian ele- 
ments n and i are 


0( a f • E 








IS;S, 


(3-57a) 


(3-57 b) 


A A (l-p^-a^)T 

t a WE - — 7 7- » C; 


(3-57c) 
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where, clearly. 
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These expressions are easily obtained from the transformation matrix in 
Appendix A of Reference 5. 

3.2. 1. 1. 1 The Function in Terms of the Direction Cosines ot, T 

It follows from the definition of the direction cosines that 


vr-*» 


(3-58) 
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By direct substitution and some algebraic manipulation, it can be shown that 
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£xp [j(s-m') ^ ] 6Xp [jlm-l«))xi] » \ y i-Ti / 


Using the relations 
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(3-60a) 
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(3-60b) 
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(3-60c) 
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it can also be shown that 


-im-s ma-s 
^\|3L ^i/SL 


- a^u.r(yF^) 


,m-s 


ni-s ^Yfi S / 

s..,i = 1 (urt yi^) 

S*rr\ -5 , / - \S-yyi 

-I/I - a (i»i) yi^) 


(3-6 la) 


(3-6 lb) 


(3-61C) 


Substituting Equations (3-59) and (3-61) into Equations (3-55), using the relation 

A 




(3-62) 


and examining the direct and retrograde cases yields (after some algebraic ma- 
nipulation) the following final expression; 


2L 


I (-1) ^ (li-lt) (t) 

(-1^3 S-m) 

(3 -63a) 

i xm-5 .-m (itm)! (1-m)! ^ ^m•Is y . „ 

^ ^ tir.): »■.)! 

(-m ^ s ^ m) 

(3-63b) 


(mSsSA) 

(3-63c) 
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similarly, the collapsed form of the function S.^ ’ takes the form 
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(O^s^m) (3-64a) 


a (uit) Uti 


(m^s^A) (3-64b) 


,m, s 


3. 2, 1, 1, 2 The Function In Terms of the Equinoctial Elements p and q 

It follows from Equations (3-57) that 


>aT 






■ ■; — 5T^ 

l^p *1 


(3-65a) 


(3-65b) 


C3-65C) 


Sulistltutlng these expressions Into Equations (3-G3) and (3-64), respectively, 
and simplifying yields the expressions 
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(- i < s < -m) (3-66a) 


rn,S J J (KmV.(J-nO! , . a aO*” v\ / /•} ««k\ 

hi' \ ^ u T ri'.U-;)" Pi-m (V) (-m^s<m) (J-66b) 


l‘ (-1^'* (p^.jlc^')’'*"' (Jf) (m< s S i) (3-66c) 
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3.2.2. 1 Choice of the Expansion Variable 

Any one of the mean, eccentric, and true longitudes can be chosen as the expan- 
sion variable. In most applications, the mean longitude Is chosen since It Is a 
simple linear function of the time which Is the natural Independent variable. The 
mean longitude is therefore well suited for generating ephemerldes or any appli- 
cation where the time history of the elements Is closely monitored. 

Use of either the true or the eccentric longitude In the disturbing function with 
time as the independent variable In the equations of motion requires that more 
complicated expressions relating the chosen longitude to the time be evaluated 
at every step. If the eccentric longitude is chosen, then Kepler's equation 

A ^ F - k sinF hcos F (3-71) 

must be evaluated. If the true longitude Is chosen, then the expression 



must be evaluated In addition to the evaluation of Kepler's equation. 

Transforming the Independent variable in the equations of motion to the desired 
longitude will remove the necessity of evaluating these expressions, but the rela- 
tion to time is no longer apparent. If the explicit time dependence Is required, 
one or both of the above expressions must be Inverted, depending on the longi- 
tude. The Iverslons of Kepler's equation requires an Iterative procedure, e.g. , 
the Newton- Raphs on method, and Is considerably more expensive than simply 
evaluating the expression. Consequently, the eccentric and true longitudes as 
independent variables are best suited for those applications where knowledge of 
the time Is not required or, at least. Is required infrequently. 
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In nciditlon to the above consiclei'atlons, the particular characteristics of the 
Fourier expansion must be considered. In Section 2.2, Fourier expansions for 
the function 



e 


jfcL 


were developed in terms of the mean, eccentric, and true longitudes. For the 
nonspherical gravitational disturbing function, n is a negative integer, 1. e, , 

n a - (Ai-ll 


For negative values of n, the expansion in the mean longitude Introduces the 
infinite series 


^ V® 


(3-73) 


* 1 s 

where the modified Hansen coefficients, Y ’ , (defined In Section 2 in Equa- 

^ o 2 

tion (2-305)) are infinite power series in either |3“ or e , both of which are func- 
tions of h and k . The expansion in the eccentric longitude also Introduces a 
similar Intinlte series. 

In contrast, for negative n , the expansion in the true longitude Inti'oduces the 
finite expression 


(i) 
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u*t) 




(3-74) 


Consequently, when choosing the expansion variable, the simple time-mean 
longitude relation and the infinite Fourier series representation must be weighed 
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against the more complicated time-true longitude relation and a closed-form 
Fourier series representation. Expansion of the nonspherlcal gravitational 
disturbing function In the eccentric longitude offers the worst of both of the other 
alternatives, thus introducing Kepler's equation as well as an Infinite Fourier 
series In the development. If the complete nonspherlcal gravitational disturbing 
function is required, the true longitude appears to be the better choice for mod- 
erate to large eccentricity satellites and the mean longitude is the better choice 
for small eccentricity satellites. 


However, in the averaged disturbing function, all short-period contributions are 
to be eliminated, regardless of the particular formulation used. In the absence 
of resonance, only the constant term of the Fourier series should survive the 
av'eraging process.^ Hence, any advantage in one expansion ov'er the others 
depends ultimately on whether the constant term in one of the expansions has a 
computational advantage over the const.uit terms in either of the other two expan- 
sions. This discussion in Section 2.2 gives simple relations between the constant 
terms in the three Fourier expansions (see Equations (2-211) and (2-212)), l,e. , 


y"' w" 


(3-75) 


and, therefore, no one constant term possesses a significant computational ad- 

var.tage over the other two. Thus, in the absence of resonance, Insofar as the 

av'eraged nonspherlcal gravitational disturbing function is concerned, it makes 

2 

no difference whether the true or mean longitude is chosen for the expansion. 


^Thls, however, requires certain assumptions about the perturbation model. 

This is discussed in more detail in Section 3.3. 

0 

“However, for the formulation of the first-order short-period variations In the 
osculating elements (Reference 5, Section 4), which requires the short-periodic 
pax*t of the disturbing function, the finite formulation afforded by the true longi- 
tude expansion Is very Important In all but the very-near circular cases. 


0 

3 

0 

0 

'! 

l). 

1 ; 

I: 

[ 

s 

i 

I 

I 

0 

0 


u 


3-26 



— / I I 1 I 1 !...' i_J_a— -.-1 


.t 


u 


Li 


S I 


. i 
K I 


^.i 




I "f 


1 ; 


jLi 


If the posstbllty of resonance phenomena Is considered, the choice of the expan- 
sion variable is no longer arbitrary. The mean longitude (anomaly) must be 
selected for the Fourier series expansion if the resonant contributions are to 
be isolated. Tills is because a resonance occurs when the ratio of the mean 
motion of the satellite, n, to the central body rotation rate, to , is very nearly 
the ratio of two small integers, l.e. , 


n 


N 

w ^ N' 


(3-76) 
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This causes the term 


N'A - N8 « a(t) 


(3-77) 
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I'eferred to as the critical term, to vary quite slowly, thus introducing a very 
long-period, A -dependent component to the motion instead of the usual short- 
period contribution produced in the absence of resonance. There is apparently 
no corresponding formulation of the critical argument in the true or eccentric 
longitudes (anomalies). 

3.2.2. 2 Introduction of the Fourier Expansion in the Mean Longitude 

Since the averaged equations of motion are to be developed for the resonant 
tesseral harmonies in addition to the zonal and uonresonant tesseral harmonics, 
the Fourier expansion in the mean longitude is required. Substituting Equa- 
tion (3-73) into Equation (3-68) yields the following complete expansion of the 
disturbing function: 

*■ 1 Z Z S Z 


i«i m«0 %«•! 


(3-78) 
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The cox'respondlng collapsed expression for the disturbing function, where 0< s<A 
and O^q < 00 , follows Immediately by replacing s by 6s and q by Vq in the above 
expression, where the parameters 6 and V are mutually Independent and assume 
the values 

( -1 (s<0) 

6 » { (3-7 9a) 


-1 

(S<0) 


(s>0) 

-1 

(q<0) 

¥l 

(q>0) 


(3-7 9b) 


and the quantities s and q are restricted to be nonnegative integers. In view of 
Equation (2-52), the resulting expression for the collapsed disturbing function Is 


m»0 s»0 0*0 6**1 V*tt 

U.»v. e J 




(3-80) 


where Sg and 6q are defined by 
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(3-81) 


The above expression Is simplified by substituting the values of 6 and V In the 
factors 
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and using the definition (Equation (2-305)) 


• (k-jnh) 


whei'e T|= sgn(Vq- €s) and 
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(3-83) 


(3-84) 


li 


which follows fi'om Equations (2-134) and (2-303). The final result can be 
expressed in terms of € and V as 


5 -i.-l.fis 


= (k-inh) 


S e 


(3-85) 


where TJ = sgn( g (V q - s)] . 


Substituting the right-hand side of Equation (3-85) into Equation (3-80) yields the 
final expression for the collapsed nonspherical gravitational disturbing function, 
i. e. , 


m»0 »»0 fi»*t v»*l 

Cits ty*nl 


(3-86) 
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where = 4 sgn(vq - s) . 
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3,3 THE AVERAGED DISTURBING FUNCTION 


This section discusses the application of the averaging operation to the nonspher- 
Ical gravitational disturbing function. The concepts of time-independent and 
time-dependent averaging are discussed and related to the classical assumptions 
of a stationary central body and of exact resonance. The discussion of resonance 
Is extended to Include the phenomenon of near resonance. 


Also presented are the time independently averaged disturbing functions for the 
zonal harmonic, the combined zonal and tesseral harmonic, and the reduced 
resonant tesseral harmonic fields. 

3.3.1 Application of the .Averaging Operation 

The flrst-oi'der averaged equations of motion require that the disturbing function, 
expressed in terms of the mean elements (a, X ) , be averaged over an appropri- 
ate Interval to remove all short-period components. The development of these 
equations was discussed In Section 3.2 of Reference 4. In the absence of reso- 
nance phenomena, the appropriate form of the averaging operation Is 


< 





^ / Rd.A) dX 
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(3-87) 


However, In the case of resonance, this averaging operation Is strictly valid 
only for properly reduced force models consisting of quasi-isolated resonant 
terms in the disturbing function.^ If the force model Is not properly reduced, 
the averaging operation should be defined as 
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R(t,A) dA 


(3-88) 


-M-n 

where N Is the number of satellite revolutions performed during the smallest 
period common to both fast variables (In this case, the satellite mean-mean 


1 


See Section 3.4 of Reference 4 for a detailed discussion. 
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longitude and the Greenwich Hour Angle or Its equivalent describing the central 
body rotation). For example, for the case of 2:1 resonance (l.e. , n/(0 - 2/1), 

N = 2 and for the case of 7 :4 resonance, N = 7 , 

In practice, the averaging Interval will not generally be centered around the origin 
of longitudes, l.e., A =0, as Implied by Equations (3-87) and (3-88), but will. In 
fact, be centered about the value of the mean-mean longitude associated with the 
time of the numerical integration step, A^. Consequently, Equation (3-87) Is more 
correctly expressed as 


Ao> 7 T 

I R(S,J)dA 

V f\ X-n J 


(3-89) 


The same modification is also appropriate for Equation (3-88). 

This dependence on the value of the mean-mean longitude at the Integration step 
time has Important implications, particularly for the numerical averaging method. 
This is discussed In Section 3. 3. 1.2. 

One form of the averaged disturbing function Is obtained by applying Equation (3-89) 
to Equation (3-78), which yields 
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clearly, only the imaginary exponential function is dependent on the mean-mean 
longitude, A, and, therefore, 
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3.3, 1. 1 Evaluation of the Averaging Integral 

The evaluation of the above definite integral is straightforward. The classical 
approach has been to assume that either Q is completely independent of A , 
which tacitly requires a constant Greenwich Hour Angle for Earth satellites, 
or that the ratio of the central body rotation rate to the mean-mean motion of 
the satellite is U; the ratio of two Integers, i.e. , exact resonance. 

The first approach of holding 9 constant over the averaging Interval during the 
averaging operation will be referred to as time-independent averaging. How- 
ever, since the value of 0 is to be evaluated at each integration step, the effects 
of the rotation of the central body are Introduced in the long-term satellite mo- 


The second classical assumption of exact resonance is a special case of time- 
dependent averaging. Specifically, the Greenwich Hour Angle, 9 , (or its equiv- 
alent for some central body other than the Earth) is permitted to vary during 
the averaging operation according to the constraint 


(3-92) 
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or, equivalently, 


]A - m0 » 0 


(3-93) 


where the integers j and m ai*e multiples of the Integers N' and N, respec- 
tively, l.e.. 
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(3- 94a) 


(3-94b) 


In practice, neither of the classical assumptions may be strictly valid for a par- 
ticular satellite. In the following discussion, the more general time -dependent 
averaging approach Is used in evaluating the definite integral in Equation (3-91). 
Each of the special cases corresponding to the classical assumptions is then 
deduced from the general result. Viewing the special cases against the back- 
ground for the general result provides additional insight into the application of 
the method of averaging. 

The method of averaging requires the disturbing function to be average.i over 
some time interval, in this application, the mean revolution period of the satel- 
lite, Because of the simple relation between the mean-mean longitude, X , and 
the time, this requirement Is easily translated to averaging over A on the Inter- 
val (0, 2rrl. Hou-ever, the time dependence of other parameters, l.e., the 
Greenwich Hour Angle (or its equivalent) should not be discounted simply because 
the averaging operation is defined in terms of A instead of time. 

mm 

The Greenwich Hour Angle (or its equivalent) is easily expressed in terms of A 

mm 

to accommodate the A-form of the averaging operation. First, 0 can be ex- 
pressed as a function of time through the relation 
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where OJ Is the rotation rate of the central body, considered to be constant ove. 
the averaging interval, and is the apparent Greenwich Ho» ■ Angle at the 
integration step time. Equation (3-93> Ignores the effects of precession and 
nutation on the value of 9 , over the averaging interval, except at the center of 
the interval. 

For the purposes of the following discussion, it is assumed that the ratio (J/n 
is constant with respect to the time. This assumption is not inconsistent with 
the basic assumption of the method of averaging, i.e. , that the slowly varying 
elements remain constant over the f veraging interval. 

Tlie mean-mean longitude of the satellite is expressed explicitly In the time by 


A * nt ■►X, 


where n is the mean-mean motion of the satellite defined bv 


(3-96) 
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The quantity jul is the gravitational parameter of the central body and a is the 
mean semimajor axis. Hence, the mean motion, n, is either constant or 
slowly varying, depending on the exact perturbations. 

Eliminating the time between Equation (3-95) and (3-96). i.e.. 
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vie Ids the relation 
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The special cases of a constant Greenwich Hour Angle and exact resonance cor- 
respond, respectively, to the conditions 
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where N and N' are Integers. 

Substituting Equation (3-99) into the definite integral in Equation (3-91) yields the 
expression 



(3-101) 


. Evaluation of the right-hand side of Equation (3-101) yields the general result 
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(for q = m = 0 or 
q - mu)/n) 




(otherwise) 
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Inspection of this I’esult Indicates that the time-dependent averaging operation 
does not generally remove all satellite-dependent (^-dependent) 8hort~period 
terms In the disturbing function. Nor does It I'emove the 0 -dependent medlum- 
and short-period contributions except in the case of resonance. (This last 
observation Is, or course, to be expected, ) The residual short-period terms 
shown in Equation (3- 102b) possess the periods 


T 

or? - mu) 


(3-103) 


where Z Is the mean period of the satellite, l.e. , 
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(3-104) 


Their contributions at the Integration step time, 1. e. , A - Aq , may survive the 
averaging operation. The sine factor In Equation (3- 102b), referred to as the 

« 

averaging factor, determines whether In fact the I'esldual short-period effects 
I'eally persist after the averaging operation and. If so, their degi'ee of signifi- 
cance. 

3, 3, 1. 2 The Averaging Factor 

It Is apparent that the averaging factor acts to suppress both the amplitudes of 
the residual sate Hite -dependent short-period terms and the amplitudes of the 
0 -dependent medium- and short-period terms, since 
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Is always satisfied, Tlie degree to which the averaging factor actually suppresses 
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the I'eslclual periodic terms depends, among other things, on the period of the 
term. 

Inspection of the right-hand sides of Equations (3-101) and (3-102) shows that 
frequency of a given term is proportional to the denominator of the averaging 
factor. Hence, the averaging factor is more effective in suppressing shorter 
period terms in general. The effect of the averaging factor is discussed sep- 
arately for the zonal and tesseral harmonic terms. 

3.3. 1.2. 1 The Averaging Factor for the Zonal Harmonic Terms 

The zonal harmonic terms in the disturbing fr.nctlon are those for which m = 0 . 
Consequently, the averaging factor in Equation (3- 102b) reduces to 


1 (for q ^ 0) (3- 106a) 

s»n(Q7r) 

i — a 0 (for q > 0) (3- 106b) 
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The case where q = 0 corresponds to those terms which contribute only to the 
long-period and secular motion of the satellite. The case q ^ 0 corresponds to 
the A -dependent short-period terms wiiich are completely suppressed by the 
averaging factor. 

3.3. 1.2.2 The Averaging Factor for the Tesseral Harmonic Terms 

The tesseral harmonic terms (Including the sectoral terms) in the disturbing 
function are those for which m ^ 0 . Tl\e iHscusslon for the tesseral harmonic 
terms is presented separately for the m-dally nonresonant terms (q - 0), the 
general nonrcsonant terms, and the resonant terms for which 


is satisfied. 
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3.3. 1.2.2. i The Averaging Factor for the m-Dally Terms 


Tlie m-dally terms In the tesseral harmonic field are those terms for which q = 0. 
Tlie periods of the m-daily terms are given by 


(3-107) 
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where t' Is the rotation period of the central body. Hence, the m-dally terms 
produce effects in the satellite motion with a frequency of m cycles per day 
(rotation period). 

Substituting q = 0 in the expression for the averaging factor in Equation (3- 102b) 
yields the averaging factor for the m-dally terms, which is 
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Clearly, 
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Equation (3-109a) is easily obtained by using the theory of limits on Equation (3-108) 
or substituting the values q = 0 and o/n - 0 directly- into the right-hand side of 
Equation (3-101) before evaluating the definite integral. Equation (3- 109b) follows 
from the pi'opertles of the sine function. 

Equation (3- 109a) corresponds to the classical assumption of a stationary central 
body. In view of Equations (3-109b), the significance of the m-dally terms in the 
averaged dlstiu'blng function is inversely proportional to the rotation rate of the 
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central body and directly proportional to the mean-mean motion of the satellite. 

It is not surprising that the m-dally effects become less significant as the central 
body rotation period approaches the satellite moan revolution period. Further- 
more, for supersynchronous satellites, the m-daily effects become much less 
significant than the -dependent short-period terms. 

In addition, for a fixed ratio uVn, Inspection of Equation (3-109) indicates that 
the m-daily terms become less significant for increasing order, 1. e.. Increasing 
m. This is expected, since the higher order m-daily terms are of shorter per- 
iods. 

For close- Earth satellites, the difference in the two averaging factors (Equa- 
tion 3-109) is negligible for the low-order m-daily effects, i. e. , m = 1,2,3. 
However, the discrepancy grows di'amatically, percentage-wise, as the order m 
Increases. Fortunately, the amplitudes of these high-order m-dally terms 
decrease rapidly, Tlius, although the large percentage errors contribute much 
smaller absolute errors, the effects of each high-order m-dally term is signifi- 
cantly corrupted by using the averaging factor for the time-independent averaging 
theory. 

More importantly, there exists a cutoff value of m where the two averaging fac- 
tors produce a discrepancy in the sign of the term, thus introducing a phase 
error of ir radians. This is easily demonstrated as follows. If it is assumed 
that the ratio w/n is bounded by the reciprocals of the integers k and k+1, l.e. , 
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\ - 1 for m > k 
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3,3. 1,2. 2.2 Nooresonaat Tesseral Harmonic Terms Excluding the m-Dallj' 
Terms 


Excluding the m-dally terms, the remaining nonresonant tesseral harmonic 
terms are those for which q 0 and m 0 . The corresponding averaging factor 
takes the general form given In Equation (3- 102b), This can also be expressed 
as 



(3-110) 


For the classical assumption o)/n = 0 , the averaging factor reduces to 


Q7T 

a 0 


Thus, It suppresses all residual nonresonant, A -dependent, short-period terms In 
the averaged disturbing function analogous to the case of the zonal harmonic terms. 


In the general case for which to/n ^0 , the averaging factor In Equation (3-110) Is 
appropriate and can be expressed In the form 




(3-112) 


3.3. 1.2. 2. 3 Resonant Tesseral Harmonic Terms 
Exact Resonance 

For the case of exact resonance, l.e. , where 
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the averaging factor given in Equation (3-110) takes the form 
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Those terms in the disturbing function for which 
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(3-113) 


(3- 114a) 
(3- 114b) 


produce the resonance phenomena. Tl’ j averaging factor for this limitirng case is 
unity . 

Tlie discussion for the averaging factor in the case of exact resonance is broad- 
ened to include not only the pure .esonant terms, but also the quasi-isolated and 
embedded resonant terms discussed in Section 3.4 of Reference 5. 

The entire nonspherical gravitational disturbing function can be considered to be 
an embedded resonant term in a general sense. The expansion of the disturb- 
ing function in the orbital elements results in a proliferation of terms as shown 
in Section 3. 2. More particularly, the imaginary exponential function which 
contains all the periodic information depends generally on two Indexes q and m . 
The index q is introduced by the expansion in the mean longitude, and the index 
m is the order of a given spherical harmonic term in the nonspherical gravita- 
tional disturbing function. These Indexes can assume any value in the range 
-oo<q<oo and where X is the degree of the given spherical har- 

monic term. 





A quasi-isolated resonant term Is an 3 ’ term or group of terms In the disturbing 
function for which Equation (3-114b) alone Is satisfied. Tlie range of the Indexes 
in the Imaginary exponential function Is then -oo < q < oo and m = kN (where 
N Is called the order of the resonance). Quasi-Isolated resonant terms are a 
subset of the embedded resonant term. The pure resonant terms are those 
quasi-isolated resonant terms for which Equation (3-114a) is also satisfied and 
are therefore a subset of the quasi-isolated resonant terms. 

For quasi-isolated resonant terms, the averaging factor assumes the form and 
values 

( 1 -tor 0 • IcM* 

s»n WnOtt ) 

’ |o icr 

This result Is obtained by substituting Equation (3-114b) Into Equation (3-113). 
Consequently, the pure resonant terms for which q = kN' are completely trans- 
parent to the averaging operation and the averaging operation completely sup- 
presses all other quasi-isolated resonant terms for which q kN*. 

For the remaining terms in the embedded resonant term, 1. e. , the disturbing 
function, the averaging factor retains the form given in Equation (3-113). It 
will tend to reduce the effects of the corresponding residual short-period terms 
but will never completely suppress these contributions. 

However, if the averaging operation defined in Equation (3-88) and centered at 
Xq is used, it is easily verified that the corresponding averaging factor takes 
the form 
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Clearly, this averaging factor vanishes for all cases where the argument of the 
sine function assumes a nonzero Integral multiple of tt , 1. e. , if 
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(for k 0) (3-117) 


u 


u 


This condition Is satisfied by all possible values for the Indexes q and m, with 
the exception of those that yield k = 0, which specify the pure resonant terms. 
Consequently, the averaging operation based on Equation (3-88) completely sup- 
px'esses all nonresonant terms In a general embedded resonant "term. " 

This fact Is of llttlo significance for analytical theories In which the pure reso- 
nance terms are easily Isolated. However, the significance cannot be over- 
stated for numerical averaging applications. Failure to pi'operly reduce the 
force model to only the quasi-isolated terms requires an Increase by a factor 
of N In the cost of the avex'aglng operation. If all residual short-period terms 
are to be completely suppressed. 

Near Resonance 

Exact resonance Is a limiting case which Is seldom, if ever, encountered. 
However, near resonance, defined by the consti'aint 
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Is of considerable practical Importance. Substituting the condition of near reso- 
nance Into the e.xpresslon for the averaging factor In Equation (3-113) yields the 
expression 
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The pure near-resonant terms are still specified by Dquations (3-114) and the 
averaging factor reduces to 
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(3-120) 


The periods of these near-resonant terms are given by the expression 
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(3-121) 


Clearly, as € appx'oaches zero, i.e. , as the resonance becomes deeper, the 
averaging factor approaches unity and the periods of the near-resonant terms 
become longer. Conversely, as £ grows larger, the resonance becomes 
shallower; the periods of the near-resonant terms grow shorter; and the aver- 
aging factor decreases in magnitude, resulting in the greater suppx'ession of 
the amplitudes of the near-resonant terms. 

The quasi-lsolateil near-resonant terms are restricted to those terms which 
satisfy the constraint given by Equation <3-ll4b) as in the case for exact reso- 
nance. Substitution of this constraint into Equation (3-11!:)) yields the averaging 
factor for the quasi-isolated near-resonant terms, which simplifies to 
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For q = UN, this exixresalon reduces to the expression given in Equation (3-120) 
for the pure near-resonant terms. 
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Generally, Equation (3-1221 Indicates that the residual short-period contribu- 
tions from the averaged quasi-isolated near-resonant term will be quite small 
since the numerator Is quite small for small € and small values of kX. For 
larger values of kX, the denominator grows larger, which helps to suppress 
the residual short-period contributions. 

In summary, the analysis of the averaging operation and, particular , the aver- 
aging factor shows that time -dependent averaging meth>o <i’» can fall to eliminate 
all short-period contributions In the averaged resonant a.id nonrescnant tesseral 
harmonic disturbing functions . In contrast, the classical assumptions eltuer of 
a stationary central body during the averaging Interval (time-independent aver- 
aging! or of an e.sact resonance do permit the complete elimination of all short- 
period terms. However, these assumptions can also produce an overestimate 
of the actual amplitudes of the long-period (resonance cases! and medium-period 
m-dally terms which survive the averaging operation. 


This circumstance presents something of a dilemma. For the analytical aver- 
aging method, the residual short -period effects can be suppressed by simply 
deleting them from the force model. However, this Is, In essence, equivalent 
to Inserting the averaging factor for the general time-dependent averaging oper- 
ation Into the time Independently averagetl tesseral hiirmonlcs models (resonant 
and nonresonant!. The method of numerical averaging Is, however, not amen- 
able to such an accommoilatlon, since the averaging factor differs for each term 
In the fully e.xpanded disturbing function ami since numerical averaging methods 
In general avoid this fully e.xpanded disturbing function representation. 

Generally, the choice for both analytical and numerical averaging methods seems 
quite clear. The nonresonant tesseral hariaonlcs terms should be deleted from 
the perturbation model, both to avoid the present dilemma and, more Importantly, 
to maximize the numerical Integration step size as discussed In Section 3.4 of 
Reference 3 or alternatively, to minimize the averaging Interval. Those contribu- 
tions can be evaluated from analttlcal formulas (Reference 49! and superimposed 
on the Integrated mean elements when necessary. 
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If, however, these averaged nonresonant tesseral harmonic terms ai’e retained 
in the perturbation model In numerical averaging applications, the error (usually 
small) in the amplitudes of the m-daily effects Introduced by the time-independent 
averaging operation is definitely preferable to the noise generated by the residual 
short-period terms ^ Introduced by the time-dependent averaging operation or to 
the- small step size required to accurately numerically Integrate these residual 
short-period contributions. 


! I 


For the case of tesseral resonance, time-dependent averaging is necessary. 
The residual short-period contributions are easily eliminated by restricting the 
force model to quasi- Isolated resonant terms alone or by averaging full force 
models over the appropriate number of satellite revolutions. 




3.3.2 Tlie Averaged Zonal Harmonic Disturbing Function 

The averaged zonal harmonic disturbing function consists of those terms in the 
averaged nonspherlcal gravitational disturbing function for which m = 0 . Re- 
sti’lcting the value of m accordingly in Equation (3-91) yields one form for the 
zonal harmonic disturbing function. If 
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for q = 0 


(3- 124a) 


for q 0 (3-2 14b) 


Sliort-period contributions which are propagated with large step sizes appear as 
spurious noise in the integrated results. 
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which asr<?es with Equ.'jtlons (3-1021 and (3-10(U. llius, the averaged zonal 
harmonic disturbing function simplifies to 
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(3-125) 


In view of the discussion in Section 2. 1.2.4, the summation over the index s for 
the range - £ £s < i is easily collapsed to a summation over the range 0 S s SX 
to give 
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where 
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The inclination function is defined by Equation 03-64) or (3-67). Evaluation of 
Equations 03-64) for m - 0 yields 
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Several simplifications can be made in Equations (3-126) and (3-128), First, 
from Equation (2-303), 
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where 
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(3-130) 


Also, It follows from Equations (2-134) and (2-303) that 
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Therefore, 
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In addition, it follows from Equations (2-303) that 
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Second, the Jacobi pohmomial representation of the inclination function for the 
zonal hvirmoaics case given in Equation (3-128) can be simplified since 
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and since it is easilv shown (Reference 26) that 
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Therefore, the Inclination function for the zonal harmonic terms takes the form 
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Substituting Equations (3-134) and (3-137) into Equation (3-126) yields the expres- 
sion 
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Since 
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(where Re(z) designates the real part of z), then 
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and Equation (3-138) Is further simplified to the form 
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since only the real part of this expression Is of Importance. Furthermore, the 
possibility of vanishing divisors 1s eliminated through Equation (3-133) and the 
following definition 
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The function Q« (i) Is a polynomial In X and Is simply the sth derivative of 
the Legendre polynomial 00 (Reference 2G), l.e., 


(3-143) 


S”bstltutlng Equations (3-133) :md (3-142) Into Equation (3-112) and defining 
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yields the following expression for the zonal harmonic disturbing function: 
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This form of the averaged zonal harmonic disturbing function will be used In the 
development of the averaged equations of motion presented In Section 3. 4. 1. 


3.3.3 The Averaged Zonal nnd Nonresonnat Tesseral Hax’monic Disturbing 
Function! 


The combined averaged zonal and nonresonant tesseral harmonic disturbing func- 
tion consists of all nonresonant terms In the averaged nonspherlcal gravitational 
disturbing function given In Equation (3-91). For the purpose of this discussion, 
it Is assumed that no resonance exists.- The disturbing function is 
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To avoid residual short-period terms In the averaged disturbing function, time 
Independent averaging Is used. Therefore, 
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(for m = 0) 
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(3-147) 
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^Tlie averaged zonal and nouresonant tesseral harmonic disturbing functions are 

combined to facilitate computational efficiency when both models are requested, 
o 

“Tlie combination of resonant and nonresonant spherical harmonic terms In the 
equations of motion Is not generally warranted because of the relative Insignifi- 
cance of the nonresonant terms and the unacceptably small step sizes they 
Impose on the numerical Integration procedure. 
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consequently, the time Independently nverhged zonnl nnd n„ 

. .. ® “oaresonant tessprni 

haimomc disturbing function takes the form 
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he meragtng factor ,s tmlty for alt nonvanlahtng terms and Is an over esttaate 
of the more reallstte averaging factor obtained from the time d a 

operation and .hlch Is given In Equation (3-XOhr n T T ‘ 
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where the sum over V has been performed and where the inclination function is 
defined by Equation (3-f)4) to be 
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Substituting this expression into Equation (3-150) yields the expression 
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since Xis Is restricted to take on only even values, the factor fe Is ahvaj'S 
unity and can be deleted from the above expression. This form of the averaged 
disturbing function will be used In the development of the averaged equations of 
motion presented In Section 3-. 4. 2. 
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3.3.4 The Averaged Resonant Tesseral Harmonic Disturbing Function 

The averaged resonant tesseral harmonic disturbing function Is Isolated from 
the averaged nonspherlcal gravitational disturbing function given In Equation (3-91) 
using the near-resouance constraint 
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In order to avoid confusion with the notation € = -1 , the symbol S is used in the 
above equations Instead of the 6 . used in Section 3,3. 1,2. 
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it then follows that 
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whei'e, most generally. 
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and Is referred to as the critical argument. 

Consequently, the resonant harmonic terms In the averaged nonspherlcal gravi- 
tational disturbing function In Equation (3-91) ai'e those for which 
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are satisfied. Thus, the averaged resonant tcsseral harmonic disturbing function 
Is given by 
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where k’ Is the Integer part of £/N. 

The argMiuent of the Imaginary exponential function is constant for exact reso- 
nance. 5'or near resommce, It assumes the role of the crltlciil argvunent In 
Equation (3-157). The averaging factor for exact resonance, which has a value 
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of unity, Is assumed above. If the resonance Is not a very deep resonance, the 
appropriate averaging factor 
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(where S is defined by Equation (3-153)) should be used. The inclination func- 
tion is defined exactly as in Equation (3-63) or (3-66), with m replaced by kN. 

In order to accommodate the user's desire to specify the particular terms of 
degree i and order m to be included in the nonspherlcal gravitational disturb- 
ing function, it is convenient to express the averaged resonant disturbing func- 
tion for a specific degfree and order as follows 
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where m satisfies Equation (3-159b). The averaging factor in Equations (3-161) 
(if used) assumes the form 
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Substltutl:!;;' the expressions for the funetlon S. ’ (Equation (3-(>:i)), the modi 
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fled Hansen eoefflelents, Vj ’ , and the appropriate averaging faetor given 
In Equation (3-lGl) into Equation (3-1(52) yields the result 
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where TJ^ - sgn(kN' -s) 
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Expressing the first summation and the first half of the second summation In 
Equation (3-164) over positive values of the index s (using Equation (2-51)) yields 
the form 
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This equation can be simplified by defining e = -1 for s<0 and € = +1 for 8>0, 
which vlelds the result 
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where iX - sgn (kN' - es) . 

This result Is Identical In form to that of the nonresoaant tesseral harmonic case 
given in Equation (3-152), with the exception of the factors 
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which are a genex’allzatlon of the equivalent factors In the nonresonant case. 

Equation (3-165) Is used in the development of the averaged equations of motion 
for the case of resonant tesseral harmonics presented In Section 3.4.3. 
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3.4 THE FIRST-ORDER AVERAGED EQUATIONS OF MOTION 

The first-order averaged equations of motion for the nonspherlcal gravitational 
perturbation are based on a modified version of Lagrange s Planetary Equations 
(discussed In Sections 2 and 3 of Reference 5). These equations expressed In 
terms of the mean equinoctial elements take the form 
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and the elements (a, h, k, p, q, X) are understood to be mean elements, n Is 
the mean-mean motion, and R is the avei'aged disturbing function. 

The disturbing functions developed in the previous sections of this report have 
been expressed in terms of the mean direction cosines of the equatorial 

z axis with respect to the equinoctial reference frame (?, g, w) , rather than in 
terms of the equinoctial elements p and q . Consequentl 3 % expressions of the 
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are required to modlf.v Equations (3-167) in order to accommodate the pax-ticular 
form of the disturbing functions. The following results are demonstrated in 
Appendix C of this document: 
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where C is detlned ns before. Substituting these expressions into Equations (3-107) 
yields the final form of the avorageil \’ariation of Parameters (\'OP) cx]uations of 
motion used in the current investigation, i.e. , 
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where n = n; A, B, and C are defined as in Equation (3-167); and 
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for any two variables x and y. 

These equations admit simplifications for certain cases. For example, it is shown 
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for the zonal and nonx’esonant tesseral (m-daily) harmonic perturbations. In 
addition, li is demonstrated that the zonal harmonics model admits the simpli- 
fication 
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The averaged equations of motion for the zonal harmonic model are presented 
in Section 3,4. 1; the averaged equations of motion for the combined zonal and 
nonresonant tesseral harmonics model are presented In Section 3.4.2; _nd the 
averaged equations of motion for the resonant tesseral harmonic model are 
presented in Section 3.4.3. The final form of the averaged equations of motion 
as given were implemented into the Research and Development (R&D) version 
of the Goddard Trajectory Determination System (GTDS). 

3.4.1 The Averaged Equations of ^lotion for the Zonal Harmonic Model 

Inspection of Equations (3-171) and (3-172) indicates that the averaged equations 
of motion require partial derivatives of the zonal harmonic disturbing function 
with respect to the elements (a, h, k, « , jS , V, A ) (the direction cosines 0(, 
j5, V are a redundant set of parameters). The averaged zonal harmonic dis- 
turbing function obtained in Section 3.3.2 (Equation (3-145)) is given b„v the real 
part of the expression 
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The averaged disturbing function is rearranged into a more optimal form for 
evaluation purposes before the required partial derivatives are obtained. 
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Ih'aclioal oonsidoiMtioas miuii't' that the summation over i be trunealeil at 
some finite value, i.e. , - s 1 s L. In aililitlon, a closer inspection of the func- 
tions in Equation (S-175) Indicates that, for s - i, 
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since (Equation (2-241))) 
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Also, for s = i - 1 , 

Thus, the index s can be restricted to the range OSs^i-2, along with the 
coniiltlon tliat 1 - s must bt' e\ en. 

It is often desirable to truncate on powers of e sini for certain cases. Since 

|(o<-j,8)* (kfj\i)* I » (esln O* (S-176) 

the range of the index s Is truncateil to 0<isSS<i-2, where S is the maximum 
power of eslnl retained,^ 

The exact order in which the t^vo summations are performei! depends on the re- 
currence relations used to evaluate the various quantities in the illsturbli\g func- 
tion. It is ilesirable that the recurrence relations be free of computational small 
lUvlsors and that they use simple starting values. Clearly, the complex polyno- 
mials (cx - j|d)^ ami (k ' jh)^ shoulil be e\ aluateil in order of ascendiitg (towers 
s in order to avoid a div ide operation whlt'h could IntriHluce a small divisor 

Silgher powers of the eccentricih’ will usually be retained since the function 
j,^-A-l,s polynomial in (1 - e-)~^' Truncating on the complex poly- 
nomial (k Jh)** is eqviivalent to truncating on the D’Alembert characteristic 
t'lisi in classical elements. 
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tor small eccentiMclties or inclinations. Similarly, the parallax factor {CLq/Ol) 
should be evaluated in ascending powers of Ji. The real polynomial Q , (i) 

possesses recurrence relations based on the recurrence relations for the asso- 
ciated Legendre polynomial, P, Jx). These recurrence relations are obtained 
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into Equations (2-232) through (2-240) (Reference 26). Tlie resulting fixed-order (s) 
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recurrence relations 
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can proceed in ascending values of Ji without any small dlu . . difficulty. The 
fLxed-degree recurrence 




must proceed in descending values of the index s in order to avoid the possibly 
small divisor 1 - . This direction is, however, not well suited for the complex 

polynomial recurrence formulas. This also applies to the recurrence relations 
for the function Kq ’ , since they are also governed by the associated 
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Leji'endro polynomial recurrence relations (In view of Equations (2-229) anti 
(2-302)) and thus are better evaluated usiiifr the recurrence on ascending \ alues 
of X . In fact, it is shown below that these functions can be expressed in terms 

of the polynomials Q. Ji) and, thus, Equations (3-178) and (3-179) directly 

X, s 

apply. 

Consequently, the fLxed-order (s) recurrence relations appear to be the best 
choice and require that the summation over t must be performed before the 
summation of s . It is easily verified that 
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The averaged disturbing function is then expressed as 
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For the purpose of software implementation, the following definitions are made: 
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which is strictly a real function. 
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Tlie real part of the avei^nged zonal harmonic dtstui^bing function (Equation (3- ISO)) 
then takes the form 
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Before the necessary partial derivatives for the averaged equations of motion 
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ai'e presented. It is necessary to elaborate further on the function ’ . 

As shown In the discussion In Section 2.2, 1.3.4, the special Hansen coefficient 
n s 

Xq ’ can be expressed as (Equation (2-223)) 

\sl 

\ ® Xo 

Furthermore, from Equation (2-229), 


-*-V ^ (i-Ol 


since only positive values of s are of Interest, the absolute value notation is 
dropped. Clearly, 
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in view of Equation (3-142). Since 
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Thus, the partial derivative with respect to h and k will be obtained through 
the chain rule, i.e. , 
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The following partial derivatives of the real part of the disturbing function given 
in Equation (3-183) are easil,v verified; 
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The summations are performed over the ranges 0 < s ^ S 1 L-2 and s+2 i i. 1 L, 
The following recurrence relations are used to evaluate the above quantities;. 

Recurrence Relations for the G and H Polvnomlals and Their Partial Derivatives 
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Recurrence Relations for the Quantity Vj,.s 
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Recurrence Relations for the Fxmction Kn ’ and dKn ’ /dx 
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was used. In this notation, the recurrence relations take the form 
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Since the averaged disturbing function is independent of A , it fellows that (L'qua 
tion (3-17 la)) 
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Also, it can now be demonstrated that 
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(Equation (3-174)) 


From Equation (3-172), 
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Multiph'ing Equations (3-187c) and (3-187b) by h and k, respectively, and taking 
the difference yields 



(3-196) 


Similarly, multiplying Equations (3-187e) and (3-187d) by oi and jd, respectively, 
and taking the difference yields 
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Multiplying Equations (3-188b) and (3«188a) by h and k, respectively, and taking 
the difference yields 
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^lultiplying Equations (3-188d) and (3-188c) by a and ^ , respectively, and 
taking the difference yields 
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Comparison of Equations (3-199) and (3-200) shows that 
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and Equations (3-174) is verilied. 

In view of Equations (3-174) and (3-194), the final form for the averaged equa- 
tions of motion for the zonal harmonic case is simplified to 
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where R Is the averaj^ed disturbing function, the elements ar« interpreted to be 
mean elements, and n designates the mean-mean motion. 

All of the requisite partial derivatives have been obtained and are given in Equa- 
tions (3-187) and (3-188). 

3.4.2 The Averaged Equations of Motion for the Combined Zonal and Xonresonant 
Tesseral Harmonic Model ^ 

The averaged disturbing function for the combined zonal and nonrcsonant tesseral 
harmonic model given In Equation (3-152) must be rearranged to maximize com- 
putational efficiency. First, the summation over 4. is truncated at L, i.e. , 

2 * X < L; the summation over m is allowed to vary through the range 
0 «m < M L; and the two summations over s are ..Mowed to vary through the 


This assumes time-independent averaging, i.e., a stationary central body dur- 
ing the averaging operation, and, thus, this model contains only the m-dally 
contribution of the tesseral harmonic model. 
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ranges 0<s<(m, s<L-2[ (where ( ] denotes the minimum value) and M+l<s 
< S < L-2 (where S < L-2 is the maximum power of e sin i retained in the expan- 
sion), The summation over s Is bounded above by L-2 for the same reasons 
presented in the discussion for the averaged equations of motion of the zonal 
harmonic model. The averaged disturbing function for the combi ,ed zonal and 
nonresonant tesseral harmonic field is given by (Equation (3-152)) 
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where [ ] denotes the minimum value. The averaging factors 1 and — 




are appropriate for the cases m = 0 and m >0, respectively. 
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Clearly, all real and complex polynomials should be computed recursively in 
order of ascending powers to avoid a divide operation and, hence, possible 
small divisors. The one exception is the evaluation of the real polynomials 
(1 - and (1 + since the order of Increasing powers is dependent 

on the sign of the factor €l. However, since the quantity ly is always non- 
negative, l.e. , 1=1 for Vs 0 and I = -1 for V< 0 , then (1 + IIT) si is always 

satisfied and there is no possiblity of a small divisor. The required computa- 

-i.-l s 

tional sequence for the functions Kq * is the same as in the zonal harmonics 
case, l.e., the recurrence uses ascending values of JL and a fixed value of s , 
The final consideration is that of the computational sequence of the Jacobi poly- 
nomials. 

Recurrence relations for the Jacobi polynomials can be found in References 26, 
47, and 48. The particular recurrence relation 
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is well suited in that only a single parameter, i.e. , the degree of the poli'nomial, 
n, varies. In addition, the starting values for the recurrence relation are quite 
simple, i.e., 
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and are easily verified from the definition of the Jacobi polynomial given in 
Equation (2-33). The requirement that the Jacobi polynomial recurrence be 
performed for increasing degree in Equation (3-205) requires that the index i 
be varied in ascending order while holding m constant in the first sum over s , 
0 < s S min(m-l, S), and that Z be varied in ascending order for fixed s in the 
second sum over s where m< s < {L-2, S] . The net result of all of these con- 
siderations requires that the averaged disturbing function be expressed with the 
summation over i as the innermost summation. 


The order of the summations over m and s depends on the range of s . The 
order of the summation'’ is rearranged as follows. First, if the following 
summations are considered: 
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it is easily -verified that 
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Thus, the net result of the above considerations on the ordering of the averaged 


disturbing function is 
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Wliile the above ordering appears to satisfy the requirements of the various 
recurrence relations, It also Imposes the recomputation of certain functions 
or, alternately, sufficient storage requirements to store these quantities. For 
the purpose of this Investigation, It Is assumed that auxiliary core storage 
resources are limited. In addition, the unnecessary recomputation of quantities 
will be avoided. 

Inspection of the second series of summations In Equation (3*-214) Indicates that 
-i-1 s 

the functions * will have to be reevaluated each time the outer sum 

over m Is Increased by one. The alternative Is to store all possible values 
for this function, which number (L-l)xS « (L-l)(L-2) . Tlie storage require- 
ment can be minimized by rearranging the outer two summations, l.e. , 
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Equations (3-214) then take the form 
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where S ^ L-2 . 
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As a result, the computation of the function ® Is not affected by the m 

summation and, hence, they are evaluated only once as they are needed. How- 
ever, the complex polynomial (a - no longer proceeds in ascending 

powers for 6l = +1, Consequently, the recurrence relation requires that they 
must be evaluated and stored prior to entering the second series of summations. 

The maximum number of quantities thus stored is 2S « 2(1-2) , which is gener- 

-Jt -1 s 

ally significantly fewer in number than the total number of polynomials, * . 

Close inspection of the summations in Equation (2-215) Indicates that 
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and the expression for the real part of the disturbing function in Equation (2-215) 
can be simplified to 
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The following partial derivatives are required for the averaged equations of 
motion due to the combined zonal and nonresonant tesseral harmonii. terms: 
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In order to obtain the partial derivatives with respect to h and k , it is necessary 
to take the partial derivative with respect to x. Thus, the complete partial de- 
rivatives with respect to h and k , respectively, are designated by 
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denotes the partial derivative with respect to h while k and x are held constant. 
In view of Equations (3-186), then 
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The pai'tinl derivatives of the disturbing function given In Equation (3-2 IT) 
require the following partial derivatives; 
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It l3 just as practical to compute the partial derivatives of the functions K, 
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and R . (t') dlrectlv via recurrence relations than to express these partial 
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derivatives in terms of the functions themselves. 
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Two linear combinations of the indexes s and m appear in braces 'in Equations 
(3-234b), (3-234c), and (3-234e). The top quantity appearing Inside the braces 
is valid for m < s and the bottom quantity is valid for m > s « In addition, the 
averaging factor Is unity for m = 0 or w/n « 0 and is 
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The following recurrence relations are used to evaluate the functions in the partial 
derivatives. 

Recurrence Relations for the Pol^vnomtals ^t. s.m 

The general recurrence relations take the form 
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These recurrence relations will be represented by the notation 
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for o.>0, b>0, and c>0. . 
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The followti i recurrence relations are used for mils 0, for which 
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The following recuri*ence I'elattons ai'e used for the case m > s >0 : 
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Reeurrenee Relations for the Functions Ej^ , \’£ ’ (Xl, 

and dx 

The structure Imposed on the expressions for the dlsturbln>j function ^Equa- 
tlon (11-2 ITR and its partial ilerivatlves (Equations (H-22dR requires that recur- 
rence relations for varying i, must be employed first, followed by the necessary 

relations for varyliyir m. The recurreiuv relations for varylUkt i used to eval- 
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to X, the proiluct \ | those proiiucts Is 

stored In a sli\jtly-iUmensloned array to provide the necessary back values re- 
quired for the recurrence relation used to evaluate the products for varyinjj 
values of m . More spcclflcalVv, the first time through Ute summations s « 0 , 
m ^ 0 , and 2 < i s L . The recurreniv relation 
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After the summation over m (0 « m ^ Mi is perfornunl. the value of 8 is increased 
b\ one, the value ot m is reset to .*ero, and the recurrence relations over i 
(maxim, s*2i < i < L) are reiveated for the new value of s and the cycle reivats, 

Kecurrence Helatlons for the IVlvnomials llV,,, ^iVi and viU^t . 

lUvurrence Relations for the pol.vnomlals n ) defined in Equations i3-222i 

are based on the recurrence relation for the Jacobi polynomial s^iven in Equation 
l3-2ivn. The recurrence is perfonued onl.v over the index I , Chancres in the 
values ot the indexes m and s are lncor|K»rated into the starting values for the 
recurrence over i . . 
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The recurrence relations for m < s are 
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The recurrence relations for m >s are 
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3.4..‘l The Averaggd Equations of Motion for the Resonnnt Tessernl Harmonics 
Model 

The averaged equations of motion for the resonant tesseral harmonics model 
were implemented In the R&D version of GTDS in such a manner as to provide 
the flexibility for the user to specify any set of spherical harmonic terms of 
degree and order (i, m) desired up to a maximum of 10 such terms. (This 
limit Is reflected In the specific dimension statements and is easily modified. ) 

No particular relationship is assumed among any of these terms. This can 
result In the recomputatlon of certain quantities common to two or more terms 
or in an increase in the cost of etmluating some of the necessary functions. 

The form of the averaged disturbing function for the resonant tesseral harmonics 
model given in Equation (3-1G5) is used for the development of the averaged equa- 
tions of motion. The expression for the averaged disturbing function for a single 
term of degree and order (1 , m) takes the form 
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The ratio A '0 = N N” Is always positive If the retrograde equinoctial refer- 
ence system Is usetl for retrograde cases. Also, only nonnegative values of m 
appear In the disturbing function. Thus, It Is easily shown that both k and kX’ 
must be positive Integers. 
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^(1) (A'Hn’)! , (l+mV. 
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In some cases, It may be desirable to truncate on the eccentricity. Therefore, 
if all powers of the eccentricity through the nth power are retained, then the 
exponents of the (k, h) polynomials must satisfy the conditions 
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in sums Rj^ and R2, and 
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in sums R^ ;mci R^ . The summation limits for the four sums are then 

[n-kN',4] 
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(3-274d) 


whore (x, y] denotes the minimum value of x and y. 

The foui sums ai*e arranged in order of increasing powers of the («, (5) poly- 
nomial. Since this exponent depends on the retrograde factor, I, the ordering 
of the sums is also dependent on the retrograde factor. For 1 = 1, the order 
chosen Is 
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and for I = -1, the order is that given in Equation (3-266), i.e.. 


R » Rq R 3^ V R 4 ) 


(3-276) 


In each case, the placement of the underlined sum in the ordering scheme is 
arbitrary because of the appearance of the complex conjugate of the (<x, poly- 
nomials which appear in the other three sums. These conjugate polynomials 
are computed independently of the other polynomials in order to avoid the possi- 
bility of a small divisor. 

Although the four sums are arranged in order of increasing powers of the (a,^) 
polynomial, inspection of Equations (3-267) indicates that the (a,|0) polimomials 
are not ordered similarly inside all summations. Specifically, for 1=1, the 
exponent of the (a,| 6 ) polynomial progresses in decreasing order in the summa- 
tion R 3 . Similarly, for the retrograde case I = -l, the exponents of the (ct«j3) 
polynomial decrease in the summation , Thus, these polynomials (the real 
and imaginary parts) are evaluated and stored prior to performing either of 
these particular summations. 

The (k, h) polynomials were not arranged in order of ascending powers because 
of a conflict with the order of the (o(, j 0 ) polynomials and because their expo- 
nents are not necessarily monotonic functions of the index s . These polynomials 
are evaluated and stored prior to the summation in which they appear. 

The partial derivatives required for the equations of motion are developed next. 
Inspection of Equations (3-267) shows that 
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The functions Iv, , ’ ai'o infinite power series In e“ or, equivalentlv, In 
h“-U~. No finite form In x = l/(l-e-l exists and, thus, It Is not Introduced 
In this model. Tlie partial derivatives take the form 
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The partial derivatives of the first sum are 
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U' the real and imaginary parts of the complex polynomial (a^/d)* and (k + jh)^ 
are denoted bv 
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and, since s-kN'>0 in Rj, Equations (2-2S1) are expressed as 
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The pjjrtlal derivatives of the other three sums are similar In form and can be 
obtained In a straightforward manner. 

Finally, the real parts of the partial derivatives are obtained as follows. If 
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whei'e I ts ai\\’ of the pai'ameters (a. t, h, k), and, similarly-, if 
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(3-288b) 


where 1^ is either a. or A , then the real part of the partial derivatives of the 
disturbing function take the general form 
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. XUj - YVj 


(3-287b) 


The recurrence relations used for evaluating the functions in the partial deriva- 
ti\es are determined next. Since this theorj- was implemented separately from 
the nonresonant theory in Section 4.2, some of the recurrence relations used 
may differ slightl^v. 


y 

y 

0 

£ 

B 

0 

0 

0 

D 

0 

B 

0 

0 

il 


3-112 


\\ 

n 



Recun’enco Relations for the (ci.,jS) Polynomials and the (U, Polynomials 

Recurrence relations for the complex polynomials (« jp) and (k *- jr^h) are 
stralghtfonvard and are special cases of the recurrence relations for the 


t,8,m 


and H polynomials illscussed In the previous section. S|)eclflcally, If 
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Recurrence Relation for the Jacobi Polvnortiial and Its DeiHvntive 


The recurrence relation for the Jacobi polynomial given In Equation (3-204) Is 

used. The recurrence relation for the derivative of the Jacobi polynomial Is 

obtained by differentiating Equation (3-204). The factors (1 + It^) and 
-Is 

(1 If) and their derivatives are computed explicitly. 

-1-1 - 1-1 Is -1-1 ts 

Evaluation of the Functions * /dh , and dK}^y» /dk 

These functions ore evaluated using the Newcomb operator method di scussed in 
Section 2.2. 1.3.3. From Equation (2-303), 
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In view of Equation (2-185), 
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Therefore, through the nth power of the eccentricity. 
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In view of Equation (2-191), this can be expressed as 
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Also, since 
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it follows that 
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(3-302) 


(It should be noted that the Integer k In the product UN’ should not be confused 
with the equinoctial element k . ) 

Tlie vipper limit of the summation over m in Equations (3-302) and (3-303) is 
increased from n/2 to (n^l)/2 to guarantee that the partial derivatives with 
respect to h and k contain all terms through order e*^. 

The Newcomb operators are computed with the recurrence relations given in 
Equations (2-196) and (2-198), They are stored In a singly-dimensioned array 
in the order In which they appear in the disturbing function. 
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SECTION' -4 - EXPLICIT THEORY FOR THK l)l8n HIUNG 
THIKD-ROPY PEU rrUDA FION 


This section presents the explicit development of the third-body illsturblns func- 
tion In terms of the equinoctial elements of the satellite. Also, the first-order 
averaged equations of m.otlon for the special case of the nonresonant near- Earth 
(Moonl satellite are presented. 

Section 4. 1 discusses the development of the thlrd-botly disturbing function from 
the equations of motion for three mutually gravitating pi'>lnt masses. Section 4.2 
describes the representation of the disturbing lunctlon In terms of ?'-c equinoctial 
elements of the satellite and the classical elements of the thU'd body. In Sections 
4. 2. 1 ami 4. 2. 2, the rotation of the satellite-dependent and third-body-dependent 
spherical harmonic functions are presented. Generally, the' formulation Intro- 
duces two Inclination functions, one for the satellite |x>sltlon and one for the 
third-body position. In Section 4.2. It Is shown that appropriate selection of 
the coordinate reference system eliminates one of these inclination functions. 

The Fourier expansloas In the mean longitude of the satellite ami the mean anom- 
aly of the third boily are intro'.Uicetl In ,he illsturblng function In Section 4.2.4. 

Section 4.:i focuses on a s|H*clal case which is applicable tcv nonresonant near- 
Earth satellites. In Section 4.o. I, the thlrd-boily tllsturblng function is ileveloped 
in terms of the equinoctial elements of the satellite and the tllrectlon cosines of 
the third boily. The relationship between this special case and the general case 
Is explored. In Section 4, .'>.2, the Fourier series expansion In the mean longi- 
tude of the satellite Is Introduced In the disturbing function. 

Section 4.4 discusses the averaged disturbing functions for the general and spe- 
cial cases dn ^‘ctlons 4.4, 1 and 4,4.2, resptvtlvelvi and the averagi**! equations 
of motion for the special case are presented in Section 4.3. In addition, the 
algorithms which were Implemented In the Uesearch and IX'velopment dlil'l 
\erslon of the Ootidard Trajectory Deterr , nation Sxsiem (G l'DS) are descrllx'd. 
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4. 1 DE\’ELOPMENT OF THE THIRD-BODY DISTOIBING FUNCTION 

Geaeralli’, a X'ariatlon of Parameters (VOP) formulation of the third-body effects 
on the motion of a satellite assumes that the distance of the disturbing third body 
from the satellite is great relative to the satellite distance from the central body.^ 
At such distances, the effects of the nonspherical gravitational field of the third 
body can be neglected and, thus, a point mass is assumed. The equations of 
motion describing the mutual gravitational attraction of three point masses in 
an arbitrary inertial reference system (see Figure 4-1) follows from the universal 
law of gravitation and Newton's Second Law of Motion. These equations of motion 
take the form 
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- 
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(4-la) 


- Of 




(4- lb) 







(4-lc) 


where r designates the position vector in an arbitrary inertial reference system 
and the subscripts c, s , and T designate the central body, satellite, and third 
body, respectively. 


^For the case of a close third body, a restricted three-body problem treatment 
(Reference 49) Is more appropriate. 
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!T is tlosirnl'lt' to trausfor tho orijiin of the reforonoo systt'm to the center of mass 
of ’.he eentriil huiK ^see I'li;ure -l-lK Hiis i.s aeeomi-Usheti thruu,s;h the translation 
specified by the etjuatlons 
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It follows that 


r-r 

& T 


The equations of motion for the satellite anil third bod.v relative to the center of 
mass of the central body are obtained by subtractlnj; Kquatlon (4-laf from Kqua- 
tions (-l-lb) and (-l-lcl. respectively, which yields 




- - |p.|») 


(4-fa) 


’v* 2 \ 

f' . - r^) 


(4 -4 b) 


Since the satellite mass, nijj , is nejt'llglble, the second term in Equation (4-4b) 
can be omitted. Thus, the motion of the disturbing third body is that of the classi- 
cal two-body problem, which Is known. ^ Consequently, this decoupling permits 


h-'quations (4-4) could have been formulated to account for the nonsphcrical grav- 
itational attraction of the central body on both the satellite and the tldrd body. 
However, this is not necessary since these effei'ts appear only in the first term 
on the right-hand side of Kquatlons (4-4) and since the purpose of this discussion 
is to obtain the disturbing function for the thlrd-lKHiy disturbing acceleration 
represented by the second term it\ Kquation (4-4a). 
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the solution of Equation (4-4a) without solving Equation (4-4b) simultaneously. 
Generally, an independently generated ephemeris for the third body is used in 
the evaluation of the disturbing acceleration in Equation (4-4a). 

Since the disturbing acceleration is the result of a conservative force, it can be 
expressed as the gradient of a potential function. It is easily verified that 


12 1 » \f\' 


(| 2 | ■ 


where the gradient operator V is defined in Equation (3-2). Thus, Equation (4-4a) 
can be expressed as 




where the disturbing function, R„ , is given by 

•J 


(|2|' [ri*) 


The first term in the disturbing potential, 1, /|S| . is called the direct term and 
reflects the direct gi'avltatlonal effects of the disturbing body on the satellite. 
The second term, referred to as the Indirect term, reflects the Indli'ect effects 
caused by the third body on the satellite through the pex'turbing effects of the 
third body on the central body. Specifically, the indirect term accounts for the 
effects on the aatelllte motion due to the motion of the central body about the 
bai'ycenter of the centra 1-body /third-body system caused by the gravitational 
attraction of the third bodv. 
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Tho expansion of Kquatton i-i-Tl Is oonslileroil next. To obtain tho expansion of 
the disturbing: fuiu'tlon In either the ^•eneral ease or the special case which Is 
considered later In this section, It Is necessary to expaml the direct part of the 
dlsturblnjt function In powers of the ratio of the saielllte ami third-body distances, 
(r' l'M.^ Tho vector etiulvalent of the L:tw of Cosines jjl' os 

|S| • 1^-r'l * Ir-r’ (4-S) 

It' the magnitude of a vector x Is tionoted by 

X « |71 (4-9) 


then, clearl.v. 




(4-10) 


since 


r • ‘ « r r* 0044* 


(4-11) 


where ^ Is the elotigatlon angle between the vectors r and r'. Thus, 


Gm-r 

A 



(4-12) 


The dlvscusslon In this section has assumed that rv r’ : however, this need not 
be the case. The following theory Is equally valid for vei'v distant satellites 
beyond the orbit of the dlsturl'lng third botiv, l.e, , r > r* , provUled that the 
expansion Is performed In powers of tr' • r) . 
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The second tactor is the well-known generating function for the Legendre poly- 
nomials (Reference 42) defined in Equation (3-24), l.e., 

p . J (f)‘ 


and, therefore, the direct part of the disturbing function takes the form 




(cos<P) 


The first few Legendre polynomials are 


* 1 


= X 


Pj^(x) a f ’I* - 7 


and so forth. Consequently, 


Gniy Gm,j. 




Gmx »* , r-f' 

— - - 6m, 
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Thus, the Pj(eos'4') term In the expansion of the illroct par- of the tilsturbins 
potential function cancels the iiullrcct part of the disturbing iwtontial function, 
ami the expansion of the disturbing potential function In powers of (r.‘ r') takes 
the form 
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(4-18) 


The object in performing the expansion Is to facilitate the development of the 
pai'tlal derivatives with respect to the satellite elements which are required 
for the Liigrange Planetary Equations. However, the point-mass term Is 
completely Independent of the satellite position. Thus, it does not contrii ute 
to the equations of motion for the satellite because the requisite partial deriv- 
atives are identically zero. Consequently, the above expression, with the 
point-mass term, Gm.j, r' , deleted, will be used In the development o‘‘ the 
disturbing functions for both the general and s{x?clal cases. 
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4.2 A GENERAL EXPANSION OF THE DISTURBING FUNCTION IN 
EQLTNOCTLAL ELEMENTS 

The expansion of the disturbing function In terms of the equinoctial elements of 
the satellite and the classical elements of the third body requires that their dis- 
tances from the central body and the cosine of the elongation angle, between 
them be expressed directly In the elements. Tlie elongation angle is considered 
first. 

It Is assumed that the positions of the satellite and third body are given respec- 
tively by the spherical coordinates (r, B, 0) and {r\ 9\ 0'), where r>0, 

0< 0< 27 t and {- tt /2) < 0<(7 t /2), etc., and are measured relative to some right- 
handed coordinate reference system with the origin located at the center of mass 
of the central body. For the purposes of discussion, the equatorial reference 
system is assumed. It follows from spherical trigonometry that 

COS0 a sm0siri0' + CO5 0 COS0' cos(0-0') (4-19) 

Application of the Addition Theorem (Reference 30) for the Legendre polynomials 
yields 

a P^(sirt0) P^(»m0') 
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4.2.1 Rotation of the Satelltte-Depondent Sphorical Harmonio Functions 

,jm0 


The surface harmonic functions P , 


must be expressed in the 


. (sin0) 

A, m ^ 

equinoctial elements. In view of the discussion presented in Section 2. 1. 2 
{Equation (2-241), the appropriate transformation from the original reference 
system to the equinoctial reference system takes the form 
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A.vn 
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(4-26) 

(lit «v«n) 

since the latitude of the satellite is identically zero in the equinoctial reference 
system. The longitude, L, is the true longitude of the satellite, measured from 
the origin of the longitudes in the equinoctial reference system. 

The inclination function is given generally by Equations (2-43), where (fl.*, i*. tJ*) 
are the Euler angles desci'iblng the orientation of the original coordinate system 
(equatorial) with respect to the equinoctial reference system. (The symbol * in 
the above expression is used to distinguish Euler angles fi'om the classical ele- 
ments H, i, ui. 

If the original cooi'dlnate system is the equatorial system, the Euler angles are 

the same as those given in the discussion of the nonspherlcal gravitational jx?r- 

terbatlon which is given in Equations i:i-53). Any of the definitions for the 

s 

inclination function S.,^* given in Equations (3-33), (3-63), or (3-66) are 
appropriate. Substituting Equation (4-26) Into Equation (4-23) yields the ex- 
pression: 




(4-27) 
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4.2,2 Rotation of the Third-Body Dependent Spherical Harmonics 

It may be desirable to express the functions C . and S . in terms of some 

m X9 m 

set of orbital elements (e.g. , classical, equinoctial, etc.). This is partlcularl3' 
true for cases involving resonance phenomena. There is no compelling reason 
to require a "nonsingular" set of elements for the third body since they are only 
parameters in the equations of motion for the satellite and do not Introduce any 
singularities. Consequently, classical elements and the associated nodal refer- 
ence frame will be used for the representation of the third body. 

The nodal reference frame is defined by the right-handed orthogonal triad 
(Xjj, yjj, Zj^) , where Xj^ points from the center of the central body to the ascend- 
ing node of the third-body orbit, A* , relative to the equatorial reference system. 
The unit vector z„ is the unit angular momentum vector of the third->-body orbi- 
tal motion, 1. e. , 

A r K r' 

* 1 


The unit vector y„ , which completes the orthogonal triad, lies in the third-body 

A 

orbital plane, 90 degrees ahead (in the direction of the third-body motion) of the 
unit vector Xjj . The relationship between the equatorial reference system 
(X, y, z) and the third-body nodal reference system (Xj^, >’j,, z,^) is shown in 
Figure 4-2. 

Applv’ing the transformation in Equation (2-24) to the surface harmonics in 
Equation (4-24), i.e. , rotating the surface harmonic functions to the nodal 
reference frame for the third body, yields 
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where o>' and f* are the argument of pericenter and the true anomaly of the 
third body, respectively. Hence, It follows that 
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( 4 - 29 ) 


Assuming the equatorial reference system as the original coordinate system, 
the Euler angles specifying the orientation of the equatorial system with respect 
to the nodal reference system are obtained from the transformation^ 


« R^(0) R^Cn*) Tg 


( 4 - 30 ) 


where r^, and r are the positions of the satellite In the nodal reference frame 
and the equatorial reference frame, respectively, and I’ and n' are the inclin- 
ation and longitude of the ascending node of the third-body orbit referred to the 
equatorial reference system. The Euler angles corresponding to the Inverse of 
the above transformation arc required and are as follows : 

fl* « 0 (4-3 la) 

I* • i' (4-31b) 

<J* * -A' (4-3 Ic) 


Clearly, the rotation 113(0) is an identity transformation. 


4-14 




4 


I 



Substituting these angles into the general expression for the function S^’^, 
given in Equations (2-43) (ivith s replaced by p), yields 


9Vl|p P-9%1 

Sjj . e 


-*n«p m.p Tn«p,*m-p 

^V/JL ^4vp 


(-A<p^-m) (4-32a) 
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^i7X 

^iVa. P^.p 

(Cv) 


(m < p < t) (4-320 


since 




S » -S 

-Jt A 


If the function is exin'essed as 

«*J» 
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»*>»p p-m ImA* »n,p 

S * I e U 

2L ^ ^ 2t 


(4-33) 
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where the definition of follows from Equation (4-32), then Equation (4-29) 

takes the form 


Gm' (i-m)! _ , . ^.v ■i'”® 


Gm' Y" 

‘ L 


X . , (4-34) 

r ±£^1 p ia) 1^“"' 

. luwii <t,p ^ ax 


pm.X 

Xtp even ) 

Comparison of this result with the result obtained by substituting Kaula's expres- 
sion (given in Equation (3-10)) into Equation (4-23) yields the relation 

F. M . (-o' P, , , (01 u;'"'' - u"’*-*' (4-35) 

^ (fi-m)i -IX ci-m)5 ax ' ’ 


where 


1-p « a(^ 


(4-36) 


This completes the rotation of the third-body-dependent spherical harmonic func- 
tions to the nodal reference system of the third body. 

Substituting the complex conjugate of Equation (4-34) and Equation (4-26) into 
Equations (4-25) yields the following expression for the disturbing function: 


00 X 


iiili i±Pl: /^\J 

U-mV, 


i*a «n«0 $»-X p»-X 

€<1*^ } Xsp ftVCA 3 


^»*»»» p-tm m 
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»P j{su-[p(w%;')+ mfl']} 


(4-37) 
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4.2,3 Eltmination of One Spherical Harmonic Rotation 

A judicious choice of the Si’Stem of reference used for the decomposition of 
the elongation angle (Equation (4-19)) in the application of the Addition Theorem 
permits the elimination of the rotation of one set of spherical harmonic func- 
tions in Equation (4-20). The selection of the equinoctial reference system as 
the system of reference eliminates the rotation of the set of satellite-dependent 
spherical harmonic functions. Selection of the nodal reference frame, associ- 
ated with the orbit of the third body, as the system of reference eliminates the 
rotation of the third-body-dependent spherical harmonic functions. In essence, 
the argument of the single remaining inclination function Is the mutual inclina- 
tion between the Instantaneous planes of the satellite and third-body orbits, 
respectively. 

Classically, the plane of the disturbing body has been adopted as the reference 
plane for those formulations involving expansions in the mutual inclination. 

Tills choice possesses the advantage of obviating the necessity of referring the 
orbital elements (position) of the third body (usually obtained from an independ- 
ently generated ephemeris) to another reference system. Because of this con- 
sideration, the plane of the third body is adopted as the reference plane in the 
following discussion. Tlie reference system, referred to as the nodal refer- 
ence system, is shown in Figure 4-2. For the purposes of this discussion, 
the dynamical system of reference is assumed to be the equatorial reference 
system and is analogous to the original reference system in the preceding de- 
velopment involving tivo spherical harmonic rotations. 

Clearly, *’>e cosine of the elongation angle is independent of the particular 
reference system used to measure the latitude and Idngttude of the satellite, 
«P,e ) , and the latitude and longitude of the third body, ( 0 ', 9 '). In the nodal 


4-17 


I 


« -ry'i 

I ■ ■ 



t 

i 

: 

I 



f I 
I I 




♦ 

I 


r > 







\ 


•i 


refei'ent'O system, the angles <p\ O ' , tleserlblng the position of the thirti body 
on a unit sphere , are 


« 0 (-i-aSa) 

0* a <jV (4-38b) 


where 6>’ and f’ are the classical argument of pericenter and true anomaly of 
the third body, respectively’. Thus, the general expression in Equation (4-13) 
reduces to 


COStji m 


cos ^ cos [o • (<uV r)| 


(4-39) 


where the satellite latitude and longitude are measured with respect to the nodal 
I’oference system. The Addition Theorem giv en by Equation (4-23) assumes the 
form 


PAcos'^') a Re 






' m«0 

^A*m even 3 


(4-40) 


Jmd 


The spherical harmonic functions ^^^(sin<^) must be rotated to the nodal 

reference fr:ime. The Euler angles required for the specification of the function 
are obtained from the rotation of the equinoctial reference frame (f, g, w) 
(defined with respect to the equatorial reference frame) into the nodal refer*, .ice 
frame, 1. e. , 




(4-41) 
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where r^, and r^j, are the position vectors of the satellite referred to tlie noilal 
anil equinoctial frames, respectively, Ai Is the mutual Inclination between the 
two orbital planes, r' Is the lons'ltiule of the ascendlni;' noile of the satellite 
orbit referred to the nodal reference system, and T Is the angular distance In 
the satellite orbital plane between the longitude of the ascendlitg notles with 
respect to the equatorial and nodal reference systems of the satellite orbit. 

The above rotation follows immediately from Inspection of Figure 4-3. Sub- 
stituting the Euler angles 


I !' 


OJ** - -T 

-Ai 

fl* « 




(4-42a) 

(4-42b) 

(4-42c) 


Into Equation (2-43) yields 


m,s ‘••w -j - mT'l 
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(-jL<s<-m) (4 -43a) 


, , (Urn'll N 

^ ('0 77-T: TT-TT i-missm) (4-43b) 
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(U\ s V 1) (4 -43c) 
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Relationship Between the Equinoctial and 
Third- Body Nodal Reference Systems 
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The expression for the disturbing function takes the form 

00 1 ^ 




K A m’ i m') C»' 


A»3l m«0 4**i. 

(Htft 


where 


Gm' (A-mM jmtw'+A*) 

j \m * 717T 't.™'"’ ® 


Substituting the explicit express '.jns for the functions and ^ Into 

Equations (4-44) yields 


Li 

00 A A 


^ IGm'V’ V V’ 

If 


Li 

•i*X m*0 4*-JL 


C JL* 4 evert ') 

0 

(Itvrt even) 




/ r- \4 j 

P, ( 0 ) P ( 0 ) - s 

Ii4 l,m \ r'/ 34 


r \1 ">,4 ][sL'w{w'+4’)] 


Thus, a comparison with Equation (4-37) shows that one summation has been 


eliminated. However, this is accompanied by an increase in the complexity of 

ni s 

evaluating the remaining S^.’ function. 


Since the equations of motion are still referred to the equatorial system, the 
quantities t*. Ai, and T, appearing In Equation (4-43), must be related to the 
equinoctial elements of the satellite and the classical elements of the third body 
(both of which are assumed to be referred to the equatorial reference system). 
The relationship between the classical elements of both the satellite and third 
body and the quantities z\ Ai, and r are easily obtained. The relationship 
in terms of the equinoctial elements can then be constructed using the results 
given In .Appendix .A of Reference 5. 
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The relatloashlp between (r', Al, Z) and (i, A, I', A') can be obtilned from 
the formulas of spherical trigonometry or as follows. In addition to the trans- 
formation given in Equation (4-41), It Is clear from the Inspection of Figure 4-3 
that 


R^(V) R^Ca'-A) R^(.i) R^(in.)x^ 


(4-47) 


Therefore, a comparison of Equations (4-41) and (4-47) yields the relation 


R^(i‘) R^(n'-n) Rj^C-O » R^C-rO RjCiA^t) (4-48) 


Clearly, 


R^(IAvt) a a R^C'c'i • Rj(lA) (4-49) 


Substituting this result Into Equation (4-4S) and postmultlplylng each side of the 
resulting equation by R~^(IA) yields 


Rj^O') R^Cn'-n) a R^C-tO 


(4-50) 


Performing the matrix multiplications In the above expression and comparing the 
I'espectlve elements of the Uvo resulting matrices yields the relations 

tOi(A'-A) a cost' cost cosAi sir) r (4-51a) 

cos i sir)fA'-A) * cosT'smr - cosAi slur' cost (4-olb) 


* MM 



i| If 


r . f . ; 


I 


tin i 5IM (Jl'-D.) = sin Ai slnT 


cosi' sin (A'-n') a COftAi cost' iinz - ilnx' cosx 


cosi cos i' cos (A'- A) + sin i' sin i » Cos Ai cost' cost i- si nr' sin « (4-51e) 


iinAlcosT' a sini' cosi - cos i' sin i cos(a'-A) 


(4-5 If) 


•sin I sinCA'-A^ * -sin Ai sin's' 


(4-5 Ig) 


sin i' sin (A'-A^ ■ smAi sinT 


(4-5 Ih) 


COsAi = cosi cosi' -t- sini sini' cos(A'-A) 


(4-511) 


If the equations of motion ave i*eferred to the nodal system also, a new set of 
equinoctial elements (a^, h^^, q^^, Aj^) and a new equinoctial reference 

A A ^ 

frame (f w^) are introduced. In view of the definition of the equinoctial 

elements in terms of the classical elements (Reference 5, Appendix A), it follows 
that the equinoctial elements of the satellite, referred to the nodal reference 
system, are given by 


■ a. 


a e slnCuj^+lTO 

« ft cosCtOj^ + IrO 


(4-52d) 





il ^ 

I / [ i 


'H "■'[ ■■■] • -I • f '• » 

., ,i . 1 ,-L J L. 1- — 


S cost' (4-52e) 

a i, XT^ (4-52f) 


where 


0)j^ a 03 - r 


(4-53) 


Is the argument of pericenter referred to the nodal reference system and to Is 
the argument of pericenter referred to the equatorial reference system. 

AAA 

The equinoctial reference frame described by the ox*thogonal triad (f g w ^ 

Is based on the orientation of the satellite orbital plane relative to the nodal sys- 
tem . (This Is In contrast to the previous case where the equations of motion 

AAA 

were referred to the equatorial system thx'ough the orthogonal triad (f, g, \v) des- 
cribing the orientation of the satellite orbit relative to the equatorial reference 
system. ) 

The rotation from the nodal reference system to the new equinoctial reference 
system Is (see Figure 4-4) 

R^C-ItO (4-54) 


where I Is the retrograde factor. This rotation Is Identical to that expressed in 
Equation (A-4) In Appendix A of Reference 5, If the symbol fl Is replaced by T’ 
and I Is replaced by Al. 


The only classical elements dependent on the system of reference are l.n, , 
which describe the orientation of the osculating ellipse with respect to the system 
of I'eference. Tlie classical elements a, e, and S. arc Invar' ant with respect to 
a change In the reference system. 
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and the Equinoctial Reference Frame Referred to the 
H Third-Body Nodal Reference System (f^, g^, w^) 
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The Eulei* angles foi* the fimctton S,„’ correspond to the inverse of the above 
ti'ansformation and are 


OJ s -T 


i * 


n* * Ir' 


Substituting these into Equation (2-49) yields an expression for the function 
which is identical in form to that given in Equation C>^“55) , with, of course, the 
conditions 


n. * r' 
i « A ( 


In addition. Equations (3-G3) and (3-G5) are also valid, provided that (« , p, <*) 
•u'c understood to be replaced by (a^, Tf^) in Equation C'i“d3) and that 
(p, q, f ) are replaced by (p^, q^, f^) in Equation (3-GG), where 
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For this cast’ In which the equations of motion are referred to the nodal refer- 


ence frame, the disturbing function Is of the same form as that given In Equa- 
tion (4-46). The S™’ ** function arguments (ot^, ) or (Pj^, qj^, X'j) are 

consistent with the integrated elements, thus avoiding the transformations in 
Equations (4-51). Since the integrated equinoctial elements are referred to 
the nodal reference sj'stem, they must be transformed if equinoctial elements 
referred to the equatorial system are desired. (However, this transformation 
is applied only at the output points rather than at every integration step. ) Equa* 
tions (4-51) provide the basis for this transformation. 

4.2,4 Introduction of the Fourier Series Ex^^ansioos 

Tlte next step in the expansion of the disturbing function Is the Introduction of 
Fourier series expansions for the products 


I 

r* e' 


for the satellite and third body, respectively. The thlx*d-botly product need not 
ntoessarlly be expanded. However, if a properly reduced force model for thlrd- 
botly resonance cases is desired, the expansion must be performed. The con- 
siderations on wl\lch the selection of the exp'nslon variable is based are similar 
to those discussed In Section 3 for the nonspherlcal gravitational theory. Al- 
though the Fourier series expansion of the px'oduct r^ Is finite In the 
eccentric longitude, it is of no significance for the averaged disturbing function 
and is not suitable for cases of resonance. Hence, the mean longitude Is chosen. 
Also, the mean anomal^v chosen is the expansion variable to facilitate the case 
of resonance. Substituting the expansions 
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(4-58) 
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where ’ is defined in Equation (2-304) and 
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(-) ^ * 
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t»-oo 


A-l,P jU' 
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(4-59) 


where 4* is the third-body mean anomaly, into Equation (4-37) yields the follow- 
ing expression for the disturbing function in Section 4, 2, 2: 


R ^ y y V f f f lifiL p ^ 

^ a.' ^ l__, L_, L_, (4>mV. U.mM 

msO ts'l. p*-X <^«-ao i^^aQ 
^ Itp etfcn) 

/ d\A -A.i,p i [<iA-t4'-pu>'.mn‘] 


(4-60) 


(a,*) ^21 ^ 


Substituting Equations (4-57) and (4-58) into Equation (4-46) yields the expres- 
sion for the single-inclination-function form of the disturbing function obmined 
in Section 4. 2. 3, i. e. , 


00 X X CO 00 

ZIZEZ SS V- 


m«0 ^••00 t*«oo 

(its «tf«n 4 Xtni c4iM) 


/ A \4 "*> .1-1, m i [ QA -tl'- wu*>' 1 

(i) s,, X, e ^ 


(4-61) 


The expressions for the disturbing functions given in Equations (4-59) and (4-60) 
may be further improved; however, they suffice for the purpose of outlining the 
development of the disturbing function* 
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4.3 EXPANSION OF THE THIRD-FODY DISTURBING FUNCTION - A SPECIAL 
CASE^ 

Expansion of the disturbing function in terms of the direction cosines of the third- 
body position vector, relative to the equinoctial reference frame, results in a 
sbriking simplification. However, the formulation does uot lend itself to obtain- 
ing reduced force models for resonant cases. The following discussion assumes 
the equatorial reference frame for the equatit us of motion (l. e. , the satellite 
equinoctial elements are referred to the equatorial reference system). 

4.3.1 Representation of the Elongation Angle in Terms of Direction Cosines 


The definition of cos ^ in terms of the dot product of the unit position vectors 
r and r' , relative to the equinoctial reference frame, of the satellite and third- 
body, respectively, are considered. Clearly, tor the satellite. 
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0 « 0 


0 - L 


Therefore, 


eos<^ ctiQ 
C06 9 sin0 


cos L 


sm L 
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and, for the third body, 


r * 


This development Is given by Cefola In Reference 10. 


1 < I 
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where and i" are the direction cosines of the third body and are defined by 


cx S COS cos 0' 

(4-65a) 

y3 » cos^' sin 

(4-65b) 

V » 8lV}0' 

(4-65c) 


The quantities 0’ and are the longitude and latitude, respectively, of the 
third body relative to the equinoctial reference frame. 

Since 

r-r' = COS'/' (4-66) 

it then follows that 

cos'/* a (xcosL-*- PsimU (4-67) 


It also follows that cos ^ can be expressed as 


cos'/' » Re e,^ ] 


(4-68) 


where 
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OC ^ ^ S COS ^ 0 


(4-69) 
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In view of liquations (-1-G2), (4-G3), and (-t-GT), Equation (4-20) takes the form 


P^Cws4') 


P^(« COS L ^ siw L) 

r <D i* 


Re 


l_, (UmY. A,m i.m 


jm CU-6') 


(•2. >n»0 
LCtm «rtn3 


(4-70) 


It follows from Equation (4-69) that 


- jvwd* 


\ci&(p‘j * \ Vl-v^ / 


(4-71) 


Thus, 


P.Ccos‘1') s Re 
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(4-72) 
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dr 
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(4-73) 


If 




(4-74) 
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then Equation (4-72) takes on the exact form given by Cefola and Broucke in 
Equation 29 of Reference 11. 

In view of Equation (4-72), the disturbing function in complex variables takes 
the form 


— 2^ L n;;3T i?) ^ 


JlsJL msO 

Cltm ewenl 


or- 


This formulation contains only two summations In contrast to Equations (4-37) 
and (4-46) which contain four and three summations, respectively. However, 
it is pointed out that the reduction in the number of summations from three to 
two requires some preprocessing of the ephemeris data to obtain the direction 
cosines of the third body (o(, /5, V) relative to the equinoctial reference frame. 
This processing, however, is not very expensive, particularly if the ephemeris 
is in the form of Cartesian position components, since 
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(4-76a) 


(4-76b) 


(4-76C) 


4.3.2 Introduction of the Fourier Series Expansion 

i iiYiL 

The Fourier series expansion for the product r 6 is required to complete 
the expansion of the disturbing function. No similar expansion is performed 
for the third body due to the specific formulation. 
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Since this formulation ultimately assumes no resonance phenomena, there is no 
compelling need to adopt the mean longitude as the expansion variable. The 
eccentric longitude provides a finite series representation which is useful for 
the representation of both the short-period and the long-period and secular 
contributions of the disturbing function. However, the eccentric longitude rep- 
I'esentation possesses no advantage over the mean longitude for the averaged 
disturbing function. Thus, for the sake of consistency with the previous develop- 
ments, the mean longitude expansion is adopted. The resulting form of the 
disturbing function is 
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(4-77) 
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4.4 THE AVERAGED DISTURBING FL'NCTION 


The averas^'^d disturbing functions for both the general ami special cases are 
developed In this section. The averaging operation Is defined In Equation (3~S9) 
and Is repeated below for convenience: 
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(4-78) 
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4.4.1 The General Case 


tlon (4-59) yields the expression 
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(4-79) 
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This exvH'esslon considers the motion of the third body during the averaging Inter- 
val. The following development Is analogous to the correspomllng development 
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for tho nonsplierloal gravitational disturbing funotion giv»»n in Section 3.3. 


Since 


X * nt A, 


(4-80a) 


l' « n't 1- JL, 


eliminating the time in tho above expressions yields the relation 


(4-80b) 


1 * — ( X - Xq) + 1q 


(4-81) 


which is analogous to Equation (3-99) in the nonspherieal case. Substituting 
this relation into the definite integral in Equation (4-79) and evaluating the result 
vields 


±f 


X^-rr 


(for q = t = 0 or q - (4-S2a) 




(otherwise) (4-82b) 


which is the third body equivalent of Equation (3-102). 

The discussion of tlie averaging factor and the reslilual short-period terms 
given in Section 3.3 is also directly applicable for the general third-body case. 
Briefly, all the short-iwrlod contributions by the terms in the disturbing func- 
tion for which q ^ 0, t = 0 in Equation (4-79) are completely eliminated in 
both time-dependently and time-lndependently averaged third-body cases. The 
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t-monthly^ terms, l.e. , q = 0 and t 0 In EqUvitlon (4-7!)), are partially sup- 
pressed in proportion to their frequency In the time-dependent case. For the 
time-independent case, the t- monthly terms are transparent to the averaging 
operation. Tills can result In an over estimate of the amplitude of each term 
and also In a phase error of 7T radians, as discussed In Section 3.3. 1.2.2. 1. 

For the case of exact resonance, the pui'e resonant terms survive the averaging 
operation, while the other quasi-isolated terms are completely suppressed. For 
cases of near resonance, the amplitudes of the "resonant terms" are reduced, 
depending on the shallowness of resonance. Also, the quasi-isolated short- 
period terms are only partially suppressed for near resonance. A more detailed 
discussion is given In Section 3.3. The analogy between the nonspherlcal gravi- 
tational disturbing function and the third-body disturbing function Is straightfor- 
ward. 

In addition, the foregoing discussion also applies to the five-sum disturbing func- 
tion given In Equation (4-61). 

4.4.2 The Averaged Disturbing Function for the Special Case 

Application of the averaging operation given In Equation (4-77) to the disturbing 
function given In Equation (4-76) yields the expression 
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(4-S3) 


^The "t-monthly" terms In the lunar model arc analogous to the m-dally terms 
contributed by the nonsphciical tesseral harmonic model. The effects contrib- 
uted by the Sun would be "t-yearly," etc. 
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This result recognizes the time dependence of the third-body direction cosines 
and the third-body distance. However, this form is apparently not well suited 
for evaluating the above integral anal 3 d:ically. 

Kaufman (Reference 50) has developed series expansions for the direction co- 
sines up through the fifth power in the ratio of the mean motions n'/n in 
order to model the time dependence. The resulting disturbing function is then 
analid:tcally averaged. However, according to the discussion in Section 3. 3 
and its analogy to the third-body disturbing function, it is obvious that allowing 
third-body motion during the averaging operation introduces short-period terms, 
except in those cases of deep resonance. Thus, it seems that the utility of such 
expansions is limited to those cases of deep resonance. However, due to the 
treatment of the third-body position in Equation (4-83), the disturbing function 
for the resonance case is, in essence, an embedded resonant term. This also 
appears to be true for Kaufman's series representation of the direction cosines. 
Hence, the averaging operation defined In Equation (3-88), centered at Aq, must 

be used. 

• 

Of course, the direction cosine formulation could be used with the numerical 
averaging method for cases Involving deep resonance, provided the averaging 
operation (Equation (3-88)) for embedded resonant terms is used. However, 
this is clearly an expensive procedure. In addition, for such applications, the 
disturbing function in Equation (4-83) is more expensive to evaluate than the 
disturbing acceleration given by the left-hand side ef Equation (4-5). ^ 
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The time-independent averaging assumption applied to Equation (4-83) requires 
that the direction cosines and third-body distance are constant over the averag- 
ing Interval. Thus, Equation (4-83) can be expressed as 
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which simplifies to 
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Implicitly, the absolute value of the averaging factor has been set to unity and 
the disturbing function contains only the long-period and secular contributions 
and the contributions from the ’'t-monthly” terms of the disturbing function given 
by Equation (4-79). Tills Is quite satisfactory for latellltes with orbital periods 
that are much smaller than the third-body orbital period. However, for satel- 
lites with longer periods, the absolute value of the true averaging factor will 
depart significantly from unity. For the "t-monthly" terms, the averaging factor 
is glvei; by 
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where t is the index in the Fourier expansion in the third-body mean anomaly 
used for the general approach. Thus, the effects of this theory, when used on 
satellites of longer periods, will be to exaggerate these "t-monthly" effects. 
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4.5 THE FIRST-ORDER AVERAGED EQUATIONS OF MOTION FOR THE 
SPECLAL CASE OF NEAR-EARTH NONRESONANT SATELLITES 

The first-order averaged equations of motion for the third-body effects are given 
by Equation (3-167), where the elements are interpreted as mean elements. 

The disturbing function is the averaged disturbing function given in Equation (4-85) 
and, thus, the equations of motion are valid for nonresonant cases only. The 
equinoctial elements p and q do not appear in the averaged disturbing function. 
The equivalent information is contained in the direction cosines of the third-body 
position vector, in view of Equations (4-7) and the definition of the vectors 

A A A 

(f, g, w) given in Appendix A of Reference 5, i.e., 

etc. Thus, the equations of motion require the partial derivatives 

dp dot dp d/3 dp dt dp 

• 

and similarly for q . This exactly parallels the development in the nonspherlcal 
gravitation theory, which uses the direction cosines of the Inertial ^ axis with 
respect to the equinoctial reference system. In fact. Equations (3-168) apply 
also to the third-body averaged equations, with the exception, of course, that 
the definition of the direction cosines differ. This result is demonstrated in 
Appendix C. Conseciuently, Equations (3-171) also apply to the thlrd-bodi’* case. 

In addition, a simplification of the equations of motion for the nonspherlcal grav- 
itational zonal harmonic terms in Equations (3-73) and (3-74) also appears in 
this formulation for the third-body perturbation. Specifically, 



where R^^ is defined in Equation (3-172). Tills simplification is demonstrated 
below. 
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Before the requisite partial derivatives for the ecjuations of motion are obtained, 

the distui’bing function is x'eordered as follows to accommodate the recuri'ence 

relations for the functions and Q* (V): 
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then the real part of the disturbing function is given by 
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partial derivatlrrs of the cHeturbing tunctloo are 
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It follows from Equations (2-230) and (2-303) that 
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For the software implementation, the definition 


A„ . = .. -- X Q 


Ml UvlV. 


ltl,m 


Cx) 


is made. Thus, 


K, 


S,m 


m . m 

(-0 A 


and 


dK 


A,m 


iv3. 


m 


de^ 


(-U 


m dAj^ 


de' 


.l,m 


(4-101) 


(4-102) 


(4-103) 


(4-104) 


(4-105) 


This treatment of the function Kq * is somewhat different from the treat- 
ment given in Section 3. It Is presented In this manner because it reflects the 


This notation is used by Cefola in Reference 11, 
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software implementation for the third-body model. For the software implemen- 
tation, the factor (-1)^^ is Included in the definition of S , i. e., 
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(for m = 0) 
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(4- 106a) 
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Because of the commonality of the special functions between this third-body model 
and the nonspherical gravitational zonal harmonics model, the recurre-oe rela- 
tions given in Equations (3-189) through (3-191) apply. The only exception is the 
recurrence relations for Cefola’s functions and their derivatives which 
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which results In the simplification of the equations of motion, is demonstrated 
below. 

It follows from Equations (4-951 that 
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However, from Equations (4-9S), It follows that 
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APPENDIX A - DERIVATION OF VON ZEIPEL'S PARTIAL 
DIFFERENTLAL EQI'ATIONI 


Tlie derivation of Von ZelpePs partial differential equation (Equation (2-192)) 
given by 
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obviously requires that expressions for the partial derivatives e &a ’ V^e and 
z *Vbz be obtained. 
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This derivation Is based In part on notes contributed by P. Cefola. 
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SubstitutlQs the relations 
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Into Equation (A-4) yields (in view of Equation (A-3)) 
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Equation (A-7) can be expressed (after some algebraic manipulation) as 
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The expression for the partial derivative z-r is now obtained. Clearly, from 
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Equation (A-2), 
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Also, from Equation (2-265) 
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Thereiore, it follows that 
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In view of the identity given in Equation (A-8), Equation (A-12) can be expressed 
(after some simplification) as 
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Inspection of Equations (A-9) and (A-13) shows that the expressions for 
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contain similar terms. Summing these expressions .vields Von Zeipel's partial 
differential equation given in Equation (A-1). 
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APPENDIX B - A JACOBI POLYXOMLAL REPRFSENTATIOX FOR THE 
FOUPIER SERIES COEFFICIENTS s 

Tlie coefficients (Equation (2- 130b)) given below are considered first: 
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Application of the linear transformations 
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Usin? the relationship between the hypergeometriu series and the Jacobi poly- 
nomial (Reference 26), l.e, , 
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it follows that Equations (B-3) can be expressed as 
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Equations (B-5) can be simplified by considering the sign of the quantity t-s. 
For Equation (B-oa), where t-s>0, 
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Similarly, Equation (B-3b) (where t~s<0) simplifies to 
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Substituting Equations (B-8) and (B-7) into Equations (B-5) yields the Jacobi poly- 
nomial representation of Equation (B-3), l.e., 
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Substituting Equations (B-9) into Equation (B-1) and simplifying the result yields 
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i>'is expression is further simplified to 
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Tills completes the derlvtitlon of the Jacobi polynomial representation for the 
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Fourier series coefficients, * , 
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The coefficients also admit a Jacobi pohmomial representation. This 

is demonstrated using the definition given in the footnote on page 2-36, i.e. , 
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Application of the linear transformation given in Equations (B-2) to the hyper- 
geometric series in Equation (B-12) 3 'ields 
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which, in turn, admits the Jacobi polynomial representation 
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where t + s<0 In Equation (B-lia) and t-^s >0 In Equation (B-14b). This result 
can be simplified to yield 
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Substituting this last result into Equation (B-12) and siniplif\'ing yields 
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which further simplifies to 
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A comparison of Equations (B-17) with Equations (B-11) shows that 
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In addition, it is also apparent that Equations (B-11) and (B-16) satisfy the con- 
ditions (see page 2-30) 
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The Jacobi polynomial representations derived above provide a new source of 
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recurrence relations for the e\ aluation of the coefficients * . 


OKlGfNAL PAGE IS 
Ob' POOR OGAU T)- 


B-9 



^>1 .. 1 1 


•f N 

>«k«» 


APPEND LX C - DERIVATION OF THE ARTIAL DERIVATIVES 

d (q, y)/d(p, q) 


The direction cosines (0(, ^ ) with respect to the equinoctial reference frame 

A A A A 

(f, g, w) of an arbitrary unit vector v' is considered, i.e. , 
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The expression for the unit vectors (f, g, w) were derived in Appendix A of 
Volume I of this report (Reference 5) and are given by 
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The unit vector v Is Independent of the equinoctial elements p and q . (For the 
applications in this report, the unit vector v is either the Inertial z axis for the 
nonspherical gravitational theory or the unit position vector of the third body, r', 
in the third-body theory, ) It follows from Equations (C-1) that 
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Taking the partial derivative with respect to p of the unit vector f yields 
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Substituting Equations (C-7). (C^S) -inri /r o\ • * 

( ), and (C-9) into Equations (C-3a), (C-3b), 

and (C 3c), respectively, and simplifying the result yields 
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